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ABSTRACT 


This report describes the activities during the fifth six month period of the Investigation 
of acoustic propagation In the atmosphere vith a realistic temperature profile. Progress has 
been achieved in tvo major directions: comparisons between the lapse model and experimental 
data taken by NASA during the second tover experiment, and development of a model propagation 
In an Inversion. Data from the second tover experiment became available near the end of 1984 
end some comparisons hove been carried out, but this work is not complete Problems with the 
temperature profiler during the experiment hove preducsd temperature profiles that are 
difficult to fit to the assumed variation of temperature vith height, but in cases vhere 
reasonable fits have been obtained agreement between the model and the experiments are close. 
The major weaknesses in the model appear to be the presence of discontinuities in some regions, 
the low <ound levels predicted near the source height, and difficulties with the argument cf the 
Henkel function being outside the allowable range. Work on the inversion model has progressed 
slowly, and the rays for that case are discussed along wlih a simple energy conservation model of 
sound level enhancement in the inversion case. 


jK TROP VCTIpN 

During the pest six months vork on grant NAG- 1-283 hes progressed In several 
directions: the lapse cese model hes become operational and is being used to make comparisons 
vith the propagation deta obtained in the second tover experiment in October, 1 984, and vith 
data collected bg Wiltshire and Shepherd using a vlnd turbine as a source; modeling of the 
Inversion case is underway and is expected to be completed in the Spring or early Summer of 
1 985. Completion of the lapse case models was carried out by two Master's students, Alex Cheng 
and Ylplng ha. Cheng [ 1 ] used the saddle point method to carry out an approximate inversion of 
the Henkel transformed sol ution as was presented i n the previous ssmi -annual report { 2] . His 
vork Is mainly concerned vith the region Inside the shadow boundary and Ms results show the 
usual interference pattern. He did, however, use the approach of Sachs and Sllbert [3] to avoid 
the singularity occuring at the shadow boundary in the saddle point approximation and tMs 
approach ellovs the sol ution to be extended 1 nto the shadow region. Me [ 4] used approxl mate 
contour integration to extend the solution deep into the shadow region and to investigate the 
contribution of the surface wave to the total solution in the lapse case. 

The total model using both of these results shows all of the characteristics of the empirical 
modal of Wiener and Ksaat [5]: a region of Increasing excess attenuation beyond the shadow 
boundary extending to a distance beyond the boundary of the same order of magnitude as the 
distance from the source to the boundary, followed by a region of nearly constant excess 
attenuation. The model often shows an overshoot in excess attenuation in the area where the 
Increasing excess attenuation region merges with the region of ’constant' excess attenuation. 
TMs overshoot it also soon in Weiner and Keest's data [5]. The modal actually shows a slowly 
Increasing excess attenuation (as r 17 *) in the ‘constant' region In the empirical model. When 
the data of the second tower experiment became available! 6], comparison shifted to those results 
and to data obtained from the wind turbine in Medicine Bov, Montana in the Spring of 1 984. 

TMs data was obtained from Shepherd [7]. Failure of a set of temperature sensors on the 
meteorological tover during the tower experiment has caused • difficulty in curve fitting the 
assumed form for the temperature profile to the actual data but some work on tMs comparison 
has been completed and more Is underway at present. 
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Work on the Inversion cose tea progressed very slowly due to the need to complete the 
lapse case and begin comparison vlth the experimental data. However, the aquations of the rays 
have been developed and geometry of this ray diagram Is understood. This proved to be an 
extremely important poi.i In the lapse model, and realization of the need to understand this 
geometry in detail will undoubtedly prove to be very helpful In this case. 

The second section of this report describes the results obtained from the lapse model to 
dote. The inversion model is described in the third section of the report. Finally, the plans for 
the final six months of the grant are discussed in the fourth section. Copies of both theses in 
which the models were developed are included in the appendix along with a copy of the Acoustical 
Society of America paper describing the plane wave solution discussed in the previous 
Semi-Annual Report [2]. 
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RESULTS FROM THE LAPSE MODEL 


The model developed for the lapse temperature gradient case has been described In the 
Fourth Sami -Annual report 1 21, but at the time of preparation of that report the modal had not 
been used sufficiently to allow any general concl'sions to be drawn. At present the model Is 
being used to analyze the date of tie second tower experiment that was cond'rUi at 
NASA-Wallopslnthe Fall of 1984 [ 6 ]. Only a small amount ofths data collected during that 
experiment hes been modeled at present snd it Is desirable to first consider a comperlswn with 
the ampl rlcal modal of Wiener and Keast[ 5] . Duri ng that r expert merits Wiener eri Keest did not 
attempt to measure either wind or temperature 5 . «hent and, thus, the comparisons are crude 
ones, but they at least allow us to look at some results of the model. 

Figure 1 presents a comparison between the empirical model of Wiener and Keest [5] end 
results calculated by ha [4] as part of his M.5. thesis. Considering that the actual data of 
Wiener and Keest shows a scatter of typically ±1 0 dB as indicated by the shadowed region, the 
modal prediction is extremely good. The "overshooting' that the model predicts at distances 
sllghlty greater then 1 0OO m also appears to be In the data of Wiener and Keest (see Figure 2 of 
[5]). The fact that a model which does not contain the effects of turbulence can produces 
"saturating" of the excess attenuation well into the shadow region is very interesting as this 
effect has often boon ascribed to turbulent scattering. Mathematically, the saturating Is the 
result of carrying out the integration for the inversion of the Henkel transformation, through an 
approximation, for a receiver in the shadow zone. The exponentially decaying effect of rays 
tangent to the boundary produce the rapid decay near the boundary, but the effect of rays with a 

complex angle of propagation (as z-><») produce this effect and no physical explanation hes been 
given at present. 

The model has also been used to Investigate the effects of surface waves under lapse 
conditions. ha[4] found that surface waves are less likly to occur under lapse conditions than 
under isothermal ones. For example, for a geometric config^rution of the source and receiver 
with a fixed surface Impedance where a surface wave will occur In the Isothermal case. It may 
or may not occur In the lapse case, and if it does not occur in the isothermal case, it will not 
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occur In the lapM com. In addition, when o surface vtve dnes occur. It Is significant mol nig in 
the region of rapid <*-cc v just be good the shadow. Poser to the source, the much stronger direct 
end t effected waves dominate, end further from the source, the exponential decay of the surface 
wave weakens it significantly below the 'saturation" valuo of the sound level in the shedov. 

More discussion of these results is given by Me ( 4] . 

Figure 2 shows the temperature profiles ooUined during run 48, record 1 2, of the second 
War experiment [6] with a source height of 1 1 .26m and a source- receiver distance of 
1 54.3m (507.75 ft.). Also shown in that figure are five attempts to fit the experimental data. 
Two of these were least-square fits, one to the entire date set, end the second to the lower four 
Cate points. The other throe fits were obtained by subjectively choosing the three parameters to 
obtain reasonable values. Both of the least-square fits gave peculiar values of these parameters. 
Using the lower four data points resulted in an extremely large value of alpha end very flat 
profile near the ground end an extremely large value of temperature at the ground. With the 
inti re set of date a very small value of alpha was obtained. 

Only the large alpha case of the loaat-aquare fits could be calculated due to problems in the 
model and the results are shown in Figure 3 and 4 for 500 Kz and 1 000 It respective! y. The 
solid line is the result of the model, the experimental date is circled, end the broken line is the 
Weiner end Keest model, which Is just the simple source level in this case. The fact that the 
broken line does not show an abrupt decrease at small heights indicates that the shadow boundary 
has not boon reached in this case. Paak levels in the interference pattern in the 500 Hz case are 
in good agreement but the location of the maxima and minima do not agree. In the 1 000 Hz case 
the model predicts peak values about 5 dB lower than the data and again the location of maxima 
and minima are not in agreement Investigation of the ray diagram temperature parameters 
seems to Indicate that It Is unlikely that this data was taken rear the shadow boundary. Similar 
comments can be made about Figures 5 to 1 0 which are calculated with the subjective fits to the 
data. 

Figures 1 1 and 1 2 are a comparison of the data of tu experiment and the model for run 
9B, record 27 in the Mire geometric configuration m the previous case. Here only tho 
least-square fit has been used and the data was clearly taken mar the shadow boundary. In the 
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500 Hz case the model over predicts the levels while In the 1 000 Hz case agreement is quite 
good. 

Attempts to model lover sources hove shown that the dip In the velues predicted by the 
model near the source height (see the Fourth Semi -annuel report [ 2] ) is e fell ure of the model 
end further work is needed. In addition It is deer that in future experiments et least one 
additional temperature sensor is needed closer to the ground to allow good curve fitting of the 
dote. 

Attempts to model the situation existing In the wind turbine data [7] have been frustrated 
by the high source location (40 to 1 20 m) which has caused the ergument of the Henkel 
functions used in the model to be outside the allowable limits of the NASA supplied program. 

This problem was discussed by Cheng 1 1 1 but no solution was given. It appears that the model 
will have to be modified to analytically allow cancellation of two terms in the model as opposed to 
simply calculating the values and allowing the cancellation to occur in the numerical results. 
This modification M not yet been undertaken, but some preliminary work has been carried out 
to identify the source of the problem. 





INVERSION PROBLEM 

Work on the inversion problem has been progressing very slovly due to completing the 
lapse model and the need to carry out comparisons vith the experimental data In that case. The 
author vtll be only teaching half time during the spring quarter end a significant amount of 
vork Is expected to be accomplished. At present the ray equations have been developed and 
interpreted. This is the necessary first step to developing the modal as it gives both 
mathematical and physical Insight into the mathematical questions that develop vMle developing 
the actual solution to the vave equation. Lack of appreciation of this important point caused 
considerable problems vith the lapse model. The ray model also yields some interesting results 
by Itself as Is described belov. 

The rays are obtai ned by use of an acoustic form of Snell 's lav 

cos e(z) cos e s 


vhere z is 'he height above the ground of an arbitrary point and s is the height of the source, t 
is the angle the ray makes vith respect to the horizontal direction, and a is the speed of sound . 
From ( 1 ) one can calculate the slope of the rays as 


dz (■ Az+B-j 1 ^ 
— ii 1 

dr L C z ♦ D J 

after suostituting the assumed temperature profile, 

AT 

T = To, + 

1 ♦ <x Z 


( 2 ) 


( 3 ) 
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Into ( 1 ) vhara 


A= <x[ !♦ «s + AT/T^ - (1 ♦ <xs)Cos 2 0 s ] (4) 

B = [ 1 ♦ <xs ♦ AT/T^ - (1 ♦ <xs) (1 ♦ AT/TJ Cos 2 e s ] (5) 

C = <x(1 ♦ <xs) Cos 2 0 S (6) 

and 

D = ( 1 ♦ AT/TJ (1 +o<s)Cos 2 e s (7) 

Equation (2) can than be integrated to obtain the rag paths. 

In the case of an inversion *T is negative and A of equation (4) mag be positive, negative 

or zero vhile the remaining functions, B, C, and D, are positive for 6 S and z corresponding to 

phgslcallg existing rags. When A Is positive Integration of (2) glelds the function F ? ( 6 s ,z ) 
defined bg 

/ (Az + B)(Cz + D) E f A (Cz+D) ^ ^ /2 

Fi( z, 6 S ) ♦ Tanh" 1 | | (8) 

A 2 A /(A C) L C (Az+B) J 


vhere 


E = A D - B C = <x AT/T*, 1 1 ♦ <xs * AT/T J(1 * «s) Cos 2 9 S (9) 
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Note that this Is very similar to the function defining the rays In the lapse case, [2]. The A > 0 


case occurs for Cos 0 S < Cos 6, where 

r 1 ♦ AT/T^ ♦txS 1 1/2 

Cos© = | 1. (10) 

L 1 ♦ o<s J 


end represents waves that are not refracted downward t!.rough a turning point and escapa the 
trapping effect of the Inversion. 

When Cos 0 S = Cos 0. A » 0 and Integration of (2) yields 

2 

F 2 (9 s ,z) = (Cz*D) 3/2 . (,n 

3C/B 

and when Cos 6 S > Cos 6 then A < 0 and Integration of (2) gives 

/((A z ♦ B)(C Z * D)1 E f -A (Cz*D)i 1/2 

Fj(6 s , z) = * Tan' 1 | 1 (12) 

A A/l-AC] l C(AZ*B)J 
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Using these functions the rays In the Inversion case can be described. As In the lapse case 
several cases occur. 

Rays propagating upward Initially are given by 


■ = Fj( © s . 2 ) - F|(e s , s) 


03) 


where 1 ■ 1,2, or 3 depending on the Initial angle of the ray on leaving the source. Rays with i 
■ 1 do not hove a turning point end go directly from the source to Infinity. The 1>2caserays 
have a turning point at infinity. The third case, 1 -3, ere rays with a turning point at a finite 
height greater then the source height. The height of this turning point Is given by 

Z tp = - B/A. (14) 

and the horizontal location by 

r tp " ” ^j( s 

Rays initially propagating downward are described by 

r = Fj(0 s , s ) - Fj( © s . 2 ) (16) 


for i ■ 1 , 2, or 3. In all three cases the rays are reflected from the ground with the reflected 
waves given by 


r = F j ( 0 S , z ) ♦ Fj( e s , s ) - 2 Fj( 9 S , 0 ). 


(17) 
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Again in the i >1 case the raga do not have a turning point, the 1 «2 case corresponds to rags 
vlth a turning point at Infinity, and the i -3 rags have a turning point at a ft idle height greater 
than the source height and given bg ( 14). 

Ang rag vlth a turning point than goes through a succession of reflections end turning 
points. The rags vhich vere travelling upvard Initially are given bg 

r = - f 3 ( e s . z ) - f 3 ( e s . s ) -2(n-l) f 3 ( e s , 0 ) (18) 


after theg have been refracted dovnverd. Here n Is the number of turning points the veve has 
passed through. After reflection at the surface these rags are described bg 


r = Fj( 9 S , z ) - F 3 ( 9 s , s ) - 2n Fj( 0 S , 0 ). (19) 


For rays joinq dovnvard i mil ally 

r . Fj( 6 S , z ) ♦ F 3 ( 6 S , s ) - 2(n*l) Fj( 6 S , 0 ) (20) 


after reflection and 


r = - f 3 ( e s . z ) * f 3 ( e s . s ) -2n f 3 ( e s , o ) (21) 


after refraction at the turning point. Again n is the number of turning points a rag has passed 
through. Figure 1 3 is a tgpical plot of the rags using this solution. 

This solution contains one unrealistic feature, there are rags that can propagate upvard to 
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vtry large heights (Infinity In the limit) before being refracted downward to effect the solution 
on the ground. Physically one would not expect such rags to exist in reality where the inversion 
would be replaced by a lapse condition at sufficient height above the ground. These rags reach 
the ground at very large distances from the source, however, having travelled an extremely long 
path length, thus one would assume that their Influence on the t 'al solution would be minimal. 

Using this solution one can obtain an approximate result for the sound level enhoncment at 
the ground by assuming energy Is conserved between two bounding rags end comparing tl«* 
Isothermal rags to the rags for the Inversion condition. The rags which are trapped in the 

inversion are defined by -0 < 0 S < 0, where © is defined by ( 1 1 ) and is a small angle 

typically less then 10°. The acoustic intensity between this limiting ray, © 5 = 0, and the 
ground can be approximated as 


I = P 2 (2tt r z(r))/p a (22) 

where the term in parentheses is the area between the ground and the limiting ray. The height of 
this ray is z(r). 

For an Isothermal atmosphere the ray is defined by 

z(r)= (Tan0) r + s (24) 

and for the inversion case by solving ( ) 4) and ( 1 2) for z es a function of radius. Setting the 
intensities equal and solving for the ratio if the emplitudes of the acoustic pressure in the 
constant temperature and Inversion cases yields 
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r 


AT/Too 



1 


1/2 


| |( <xr ) ♦ <xs 

1 1 ♦ ocs ♦ AT/T^J 



f 3 cx r 

| /( AT/TJ ♦(! 

L 2 


ocs ♦ AT/Too ) 372 


2/3 

- I * AT/Too 


J 


( 24 ) 


and the enhancement In sound pressure level due to the Inversion relative to the Isothermal case 
Is given bg 10 Log(P|^/ P c *). This function Is plotted in Figures 14 to 1 7 for some typical 

parameters. In the limit of large radius (24) behaves like r 1 and the sound level 
enhancement like ( 1 0/3) Log (r) or about 1 dB/doubling of distance, thus the sound pressure 
level of a source in an inversion vould be expected to decrease 5 dB/doubling of distance at large 
distances from the source based on this model . It should be realized that the ground has been 
assumed to be perfectly reflecting here and an absorbing ground vould Increase the decay rate. 

Experimental data to compare to this simple model is scarce. Shepherd [7] reports a 6 
dB/doubling of distance. Wiener and Keest [5] shove positive excess attenuation in the 
dovnvlnd case, and thus a faster decay than the geometrical 6 dB/doubling of distance, thus the 
model is interesting but unsubstantiated. 
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FUTURE WORK 

During the remaining six months of the current grant, vorlc is planned for three mein 
tasks: continued comparison betveen the existing lapse model and experimental data, this may 
require some modification in the model ae discussed earlier, due to limits on the argument of the 
Henkel functions and to eliminate some of the discontinuities present in the current formulation; 
an attempt to use the fast fourier transform technique to invert the Henkel transformed 
sol ution, this vill eli mi note the disconti nui ties present i n the present r and may be 
necessary in the inversion case vhere many rays can reach an observer location, zy* 
development of the inversion model, this model is partly worked out but has not been completed 
and has not been presented here. A paper is being prepared describing the ray sol ution to both 
the lapse end Inversion problems. It is anticipated that this paper will form a foundation to 
develop both the lapse and inversion models In future papers. The author will be presaging a 
paper on this work at the Second Workshop on Long Range Acoustic Propagation in Februet y, 

1 985, and has submitted an abstract of a paper to the Acoustical Society of America for 
presentation at the Spring meeting in 1985. 
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SOUND PRESSURE LEVEL VERSUS HORIZONTAL DISTANCE 
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Figure 1. Sound pressure level versus horizontal distance for AT/T 00 = 

.001, « = 3.5 rrf 1 , o> = 10,000 sec" 1 , s = 3 m, and z = 1 m. The empirical 
model is due to Wiener and Keast [51 with <)> c = 180 *, and <f> = 0* to 

simulate temperature gradient effects. 
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Temperature - *C 

Figure 2. Temperature versus height from the second NASA tower 

experiment [G] record 12. Curve 1 is a least-square fit to the data at the 
four lowest heights, T* = 19.4 \ AT = 6.6\ <x = 21.1 m' 1 . Curve 2 is a 

least-square fit to all data except that at 1.7 m, T*, = 17.7 \ AT = 2.6*, « 

= .1 m" 1 . Curves 3 to 5 are subjective fits to the data set: T^ = 19.8 \ AT 

= 1.5 \ « = 2. m' 1 ; T*, = 19.8 \ AT = 1.5 \<x= 1. m‘ ! ; T^ = 19.8 \ AT = 

1.0 *, « = 2. m’ 1 . 
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Figure 3. Sound Pressure veiSuo height, record 12, s = 1 1.26 m, r = 154.3 
m, and f = 500 Hz. Model results with = 19.4 *, AT = 6.6’, and <x = 21.1 
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Figure 4. Sound Pressure versus height, record 12, s ■ 1 1.26 m, r - 154.3 
m, and f - 1000 Hz. Model results with T*, = 19.4 \ AT = 6.6*, and « = 


[103. Ill] 




13A31 3HnSS3«d ONflOS 


I ' 



[103. Ill] 


O 

CD 


O 

in 


o 

TT 


O 

cn 


o 

cvi 


r 

o 


13A31 3dnSS3dd ONHOS 



(103. Ill] 



"I3A33 3dflSS3Hd ONOOS 


1 


Si 


REC12 1000HZ CASE 2 


I 



13A31 3mSS3dd ONHOS 


.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 

ELEVATION 

Figure 8. Sound Pressure versus height, record 12, s = 1 1.26 m, r = 154.3 
m, and f = 1000 Hz. Model results with = 19.8 V AT = 1.5 V and « = 1 
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ABSTRACT 


A mathematical model for the propagation of acoustic waves from 
a point source in a thermally stratified medium has been developed. 
The temperature profile chosen is a model of a fully developed lapse 
condition and fits the limited amount of data available. This profile 
has a large vertical temperature gradient near the ground and asympto- 
tically approaches a constant value high above the ground. C. I. Ches- 
sers model for grass has been used to characterize the ground surface 
impedance. 

The governing equation is a simple wave equation with a variable 
sound speed and a source term. A semi -analyti cal solution to the 
problem is obtained by using a Hankel transformation and the Saddle Point 
method to invert the transformation. Some of the steps are carried out 
numerical ly. 

In the illuminated' region where geometrical rays can reach, 
an interference pattern occurs. And in the shadow region where no rays 
can enter according to the geometric acoustic, an exponentically decay- 
ing behavior is observed. The sound pressure level predicted in the 
present model is in good agreement with the empirical model formulated 
by Wiener and Keast. 
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NOMENCLATURE 


a Sound speed 

hj, h2 Modified Hankel functions 

(1) (2) 

H1/3, H1/3 Hankel functions of 1/3 order 

1 /TT 


Jo 

p 

q 

R 

r 

5 
s 
T 
t 
V 
W 

z 

z 

a 

6 
6 

AT 

0 

P 


Bessel function of first kind 
Pressure 

Source strength constant 
Reflection coefficient 
Receiver radius 
Specific entropy 
Source height 
Temperature 
Time 

Velocity 

Component of velocity along z axis 
Normal impedance 
Receiver height 

Temperature variation parameter 

Hankel transform variable 

Dirac delta function 

Temperature change between z = 0 and z 

Incidence angle 

Density 


Arbitrary independent variable 
Circular frequency 


S 

b) 

Subscripts 

M Mean values existing without an acoustic disturbance 

o Evaluated at z = 0 

» Evaluated at z » 

Supej^sc ri pt s_ 

' Variables with acoustic disturbation 



CHAPTER I 


INTRODUCTION 

The study of sound propagation has a long and intriguing history, 
and is of increased importance due to the needs of the modern society. 
With the increased population density and continuing expansion of air 
traffic, jet aircraft noise has become one of the major issues in envi- 
ronment noise control. Its impact on the airport community has led to 
the strict noise certification requirements imposed by the Federal 
Aviation Agency (FAA), which is the major impetus for acquiring a better 
understanding of outdoor sound propagation. 

A great deal of progress has been made towards understanding the 
mech<. 1 isms associated with the propagation of sound in the atmosphere. 
Four major factors are found to be significant and have been the major 
areas in recent research. They are the atmospheric absorption, reflec- 
tion from the ground surface, scattering by turbulence, and refraction 
by wind and temperature gradients in the atmosphere. Each of these 
effects is described briefly in the following except the last case, 
refraction due to temperature i nhcmogeneity in the atmosphere, which 
is the subject of this discussion and will be examined thoroughly. 

The absorption of acoustic energy during propagation through the 
atmosphere has been investigated for many years. Three mechanisms have 
proven to be the major contributors. The first one, classical absorp- 
tion, is caused by the transport process of classical physics (shear 
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viscosity, thermal conductivity, mass diffusion, and thermal diffusion). 
The second is the absorption caused by rotational relaxation of the 
molecules in air [1]. The third mechanism is due to vibrational 
relaxation of oxygen and nitrogen molecules. It was the discovery 
of the latter factor that produced agreement between the predicted 
absorption and measured data [2]. Although further investigation at 
low frequencies and lower humidity is desirable, it appears that the 
atmospheric absorption is a considerably wel 1 -understood phenomenon. 

The reflection of waves at a ground surface modeled by a finite, 
nonzero normal impedance surface has proven to be a complicated and 
intriguing subject. Several mode 1 s of the ground surface and differ- 
ent mathematical methods and approximations have been presented in a 
pletlora of papers on this subject. The 'local reacting 1 model of the 
surface is the most common and widely used in current practice, and 
in general the solution to the wave equation consists of four terms. 
The first term represents the wave travelling directly from the source 
to the receiver without contacting the ground. The second term repre- 
sents the reflected wave which appears to originate at an image source 
below the ground. The third term is generally referred to as a ground 
wave because it produces a strong signal near the ground as in the 
propagation of electromagnetic waves above the earth. This term 
becomes dominant when the receiver is far from the source and near the 
ground surface, where the direct and reflected wave terms cancel. The 
fourth term, whose exact physical nature is still obscure, represents 
a 'surface wave' [3]. It originates from the mathematics in the wave 
equation solution, and behaves as a disturbance in the air propagating 
cyl i ndri cal ly outward and decaying exponentially vith height. This 
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disturbance only occurs under some combinations of source-receiver 
geometry and values of the surface impedance. This modeling of the 
ground reflection is found to have a good agreement with the measure- 
ments conducted by Lawhead and Ruduick [4], provided that an appropriate 
normal impedence of the ground surface is available. Chessel 's model 
of the ground [5], which allows the impedance to be calculated from a 
flow resistance, coupled with the type of model described above gives 
reasonable predictions of the reflection phenomenon. 

Turbulence scattering, perhaps the least understood mechanism for 
outdoor sound propagation, is believed to be responsible for the sharp 
spikes in the measurements of an acoustic signal which has propagated 
a long distance. It is due to tne wind velocity and temperature fluc- 
tuations which cause secondary sound waves to spread out from the 
original path of the sound wave in a variety of directions. Several 
theoretical models pertaining to this subject have been considered 
over the past century. However, none of vhese theories totally explains 
the mechanisms associated with turbulence that result in the attenuation 
of sound and that can be used for a reliable prediction. 

Refraction is the bending of sound rays caused by gradual changes 
in the propagation speed of a sound wave. The variation of sound speed 
in the atmosphere is brought about principally by vertical gradients of 
temperature and wind velocity. Since the sound speed is proportional 
to the square root of the temperature, the vertical variation in tem- 
perature will cause a change in sound speed with height so that the 
wave front is moving at different speeds for various heights. This 
causes the wave to turn and is known as refraction. Similarly, wind 
gradients cause portions of the wave front to move faster as seen Dy 
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vectorially adding the wind speed and sound wave speed. The resultant 
propagation speed is not constant with neight, so again refraction 
occurs. There are two types of refractions, upward and downward. The 
upward refraction is caused by the temperature lapse condition when 
temperature decreases with increasing height, or when the sound propa- 
gates in the opposite direction of wind. Conversely, the downward 
refraction is due to the temper* cure inversion, when temperature in- 
creases with height, or when the sound propagates as the same direction 
of wind. The major acoustic effect in refraction is the production of 
a refractive shadow zone, i;he p, according to the ray picture (geometric 
acoustics) no sound may '*nter. The detailed measurements of Parkin 
and Scholes [6] concluded t .at refraction due to vertical gradients of 
wind and temperature in practice produces equivalent acoustic effects. 
However, these results are based on a relatively little information on 
the temperature condition: it is being described as lapse, neutral or 
inversion only, and thus provides no insight on the quantity of each 
individual factor. It is the purpose of this paper to investigate 
primarily on the subject of the temperature inhomogeneity in the atmos- 
phere. 

The atmosphere is not isothermal under normal conditions. During 
the daytime hours, radiation from the sun heats up the air in contact 
with the ground by means of conduction, and the atmosphere is typically 
warmer near the ground surface and a lapse condition exists. As the 
sun sets, the surface of the earth cools by radiation. This cools the 
air adjacent to the earth's surface and an inversion temperature 
gradient results. The former condition, a lapse condition with tempera- 
ture decreases with height, is of primary interest in the present 
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work. Although these regions of strong temperature gradients extend, 
at maximum, a few meters above the ground, many acoustic sources and 
receivers are also in this range of elevation. 

The influence of temperature homogeneities on sound propagation 
has been the subject of numerous theoretical studies. Quantitive model- 
ing of refraction caused by sound speed gradients has been carried out 
over the past several years, particularly in connection with underwater 
acoustics [7], but little work has been done directly related to the 
atmospheric problem. There are two closely related approaches that 
have been used, ray-tracing and a direct solution of the wave equation. 
Ray tracing provides a qualitative picture of sound behavior in the 
present of temperature gradient, but determination of the actual sound 
levels is considerably more complicated, particularly in the regions 
usually referred to as caustics or convergence zones, where the sound 
is being focused. Also, this method gives no information about the 
behavior in the shadow region. Direct solution of the wave equation 
with a variable sound speed has been undertaken by Pekaris [8], Prid- 
more-Brown and Ingard [9], Tolstoy [10] and others mainly with the aim 
of determining the behavior in the shadow region. 

A model is being developed by Van Moorhem [11] to apply to the 
atmospheric case. This work is different from the previous ones in 
that a physically realistic temperature profile and a finite impedance 
boundary condition are used to describe the ground surface. The pro- 
pagation of plane waves has already been considered by Landheim [12]. 
The point source problem, which is considerably more complex, but is 
closer to a realistic application, is considered here. Physically, 
determining the point source solution is nothing more than the summing 
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of the plane wave solutions over all possible wave incidence angles. 
However, this summing process, or integration, is expressed as an 
inverse Hankel 'transformation which poses a complicated mathematical 
problem; and no exact solution can be obtained with the existing know- 
ledge. It is at this point this investigation began. A direct numeri- 
cal integration could be carried out, but it would require a great 
deal of computation and is therefore not a practical low cost prediction 
scheme. Another alternative is by employing approximation methods; 
it is essential and indeed the only approach that can lead to an semi- 
analytic solution to our problem. This approach has the advantage 
over the numerical evaluation in that it is more useful for understand- 
ing the solution physically, although perhaps not for accuracy. 

The entire development of techniques for evaluating the formal 
solution is based on an approximation schemes in which the Hankel func- 
tions are being replaced by their asymptotic expansions. By doing so, 
the original problem is put into a form so that the Saddle Point method 
[13] becomes applicable, although some modifications are necessary. 
This approach has been used successfully in the isothermal problem and 
in underwater acoustics and would appear useful here. This method 
will be described in more detail in the following chapters. 

The purpose of this thesis, thus, is to numerically imple- 
ment this model with a computer program. With the required informa- 
tion given, such as the strength and the frequency of the sound source, 
this program is to predict the sound level received in any particular 
location in the atmosphere. The work is begun with a genera 1 descrip- 
tion of how the temperature inhomogeneity affects sound propagation. 
This is done with the ray tracing method and will be described in 
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Chapter II. Chapter III discusses the governing equation and its 
formal solution as obtained by Van Moorhem [11]. Chapter IV describes 
the analytical approach to the inverse Hankel transformation using the 
approximation (Saddle Point) method. The results and conclusions are 
discussed in the following chapters. Plots with different parameters 
are generated to give an overall depiction of the relative sound magni- 
tude versus height in the physical space and also the results are 
compared with the empirical model performed by Weiner and Keast. 
A complete list of the program and its usage is given in the appendix. 


CHAPTER II 


RAY TRACING 

Tne direction of wave propagation can be represented by rays 
which are lines perpendicular to the wave fronts. Despite the fact 
that this ray representation only provides a qualitative picture of 
the sound behavior in an inhomogenous medium, and the wave nature 
of sound is completely excluded, it explains some of the physical 
characteristics of the results obtained from the wave equation in the 
later chapters. 

The location of rays emanating from a source at a height s above 
a plane surface in a thermally stratified medium can be obtained 
through the Snell's Law [14], The vertical temperature profile em- 
ployed here is 

AT 

T ( 2 ) = Too + (1) 

1 + az 

which allows a large gradient near the ground (z = 0), but asympototi- 
cally approaches a constant high above the ground. In the present case, 
AT is a positive value which describes a lapse condition that generally 
exists during the daylight hours. For a perfect gas, the speed of sound 


is given by 

a = Art . 


( 2 ) 
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Substituting (2) into (1) yields 
o 2 AT 1 

a 2 (z) * a. (1 ♦ — "*77 — ) • (3) 

T. 1 ♦ aZ 

Snell's law requires that 

cos 0 (z) , , 

= constant (4) 

a(z) 

where e ( z ) is the angle between the ray and the horizontal in the 
present case. The constant can be evaluated by letting e(s) be the 
initial angle of a sound ray propagating out at the source location, 
and a ( s ) is the sound speed at the source height. Notice that the con- 
stant varies according to the e(s) chosen, and tnus represents rays 
with different initial angles. So, for a ray with an initial angle 
0 (s ) , the corresponding angle e(z) with the horizontal at a height 
z is 


cos e(z) = 


— — cos e (s) 
a ( s ) 


(5) 


With this relationship, between the local and the initial angle of 
incident, we can then trace the slopes of a particular ray at different 
locations. The slope of a ray is 


dz 

— = tan e (z) , (6) 

dr 


and by expressing the tangent in term of cosine 
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tan e(z) = ± / — 2 — TT _1 
'cos 0 (z ) 


we obtain 


(7) 


dz 

dr 


a 2 (s) 1 

~~T — 9 

a (z) cos 0 ( s ) 


- 1 


( 8 ) 


The choice of the sign depends on whether the ray is propagating 
upward or downward as it leaves the source. Substituting the expres- 
sions for a ( z ) in (3) and rearranging the terms yield 


dz 

dr 


AZ + B 


CZ + D 


(9) 


where 

A = a (1 + aS + AT/j - (1 + aS) COS 2 0(s)) 

B = 1 + as + AT/ t - (1 + as) (1 + AT/ t ) cos 2 9(s) 

9 do) 

C = a (1 + as) cos 2 e(s) 

D = (1 + AT/j) (1 + as) cos 2 0(s) . 

The next task is to evaluate the horizontal distance as a function of 

height. This is done by integrating (9) which gives rise to four dif- 

ferent types of rays with four correspondi ng equations given by: 

(1) r = F ( z , 9 ( s ) ) - F(s,0(s)) (11) 

for rays propagating upward from the source. 

(2) r = F(s,0(s) ) - F(z,0(s)) (12) 

for rays propagating downward from the source. 

(3) r = F(z,9(s)) + F(s,0(s))- 2F(0.e(s)) (13) 
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for rays propagating downward initially and reflected upward at the 

\ 

surface . 

(4) r = F(z,e(s)) + F(s, e(s ) ) (14) 

for rays propagating downward initially and refracted upward without 
reaching the surface. The function F(z,6(z)) resulting from the inte- 
gration of (9) is given by 


F(z,e(s) ) 


7(AZ + B)(CZ + D) 

+ 

A 


AD - BC ^ 1 + 4> ^ 

2A vfiC It - J 


(15) 


where 

C (AZ + B) 

* ~ A (CZ + D) 

These four types of rays are shown in Figure 1. With this ray tracing 
technique, it can be seen from the figure that the physical space is 
divided into seven different regions. Except for region 7, each of 
them is characterized by a combination of two types of rays. As is 
shown in Figure 1, regions 1 and 3 are both characterized by the direct 
and reflected waves, except that the former has the direct wave tra- 
velling upward while the latter is travelling downward. Region 5 is 
made up of reflected and refracted waves. It is to be stressed at 
this point that, at least from a mathematical stand point, there exist 
two kinds of refracted waves. Care must be taken to distinguish them 
because, as we will see in the later section, their mathematical 
solutions for the sound pressure are different by a 90° phase anglp. 
Physically, they are determined by whether the refracted wave has or 
has not reached the caustic before it arrives at the receiver's loca- 
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tion. The refracted wave In region 5 belongs to the latter kind. 
Regions 2 and 4 are identical to 1 and 3 except that the reflected 
wave has been replaced by a refracted wave: the one that has reached 
the caustic before it entered the region. Region 6, which is bounded 
by the caustic, consists of both kinds of refracted waves. It follows 
that at any locations, except that in the shadow region, the receiver 
is always reached by two different sound rays identified by their 
initial angle 9(s) at the source. This phenomenon, as we will see in the 
the following sections, contributes a great deal of insight to the 
physical interpretation of the mathematical solution of sound propaga- 
tion from a point source. 


CHAPTER III 


FORMAL SOLUTION OF THE POINT SOURCE PROBLEM 

This chapter summarizes the formal solution as obtained by Van 
Moorhem [11]. Sound propagation can usual’y be regarded as a small 
acoustic disturbance to the atmosphere's ambient state. Even in a 
thermally stratified atmosphere, its physical behavior is still govern- 
ed by the linear acoustic wave equations but with a variable sound 
speed. Assuming there are no viscous or body forces exerted on the 
medium, the continuity and momentum equations become 

— + V-(pV) = Q (16) 

3t 

where Q is the source term and 
3 V 1 

— + VVV = - VP (17) 

3t p 

where p, V, and P are the physical variables: density, velocity and 
pressure, respectively. 

Since the ambient state is thermally i nhomogenous , the equations 
of state for the atmosphere become [15] 

P = P (p,S) (18) 


and 
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DS 

dT 


= o 


(19) 


where S is the specific entropy of the air particle. 

For a small acoustic disturbance, the physical variables can be 
expressed as the sum of their time mean and their fluctuating parts. 
Thus assuming 

P = Pm + P'(r,z,t) 
p = pm(z) + p ' (r,z,t) 

( 20 ) 

V = V 1 (r,z,t) 

S = S M (z) 

where the quantities with subscript M represent the time mean atmos- 
pheric condition and are functions of only height above the ground 
with the exception of pressure which is assumed constant throughout 
the atmosphere. Also, the medium is assumed to be at rest (Vm = 0). 
P', p 1 , and V' represent the acoustic penturbation and always assumed 
small compared to their corresponding time mean or ambient values. Sub- 
stituting Equation (20) into Equations (16) through (19) and neglecting 
the products of the second order terms, the governing equations become 


3 P 

— + pmVV 1 + V'VpM = Q 
3t 


3 V 1 

— + — Vp' = 0 
3t pm 


3 P 1 2 

— + V'VPm = a 
3 1 M 


3 P ' 

( + V'Vp 

3t 


( 21 ) 

( 22 ) 


( 23 ) 


where 


2 3P. 

a M ' V s . s 


M 


( 24 ) 


is the variable sound speed which depends on the height z. Combining 
Equations (21) through (23) and eliminating the acoustic variables, 
with the exception of pressure, yields 


1 

M 


3 2 P 1 

at 2 


- v^p 1 - 


7p 0 VP 1 


= Q 


(25) 


Hereafter, the primes on P and V are deleted for simplicity. The first 
two terms in Equation (25) are of order (ui/aj 2 and tne third term, 
based on the temperature profile given by (1), is of order (cWa*). Here 
uj is the angular frequency of the sound source and a describes the tem- 
perature gradient at ground. Note that the term (co/aoo) is actually the 
reciprocal of wavelength and the third term is nondimensioni al i zed by a. 
Thus, if the frequency is sufficiently high that the wave length is 
small compared to the scale of the variation a in temperature in the 
atmosphere, so that w/ta^a) >> 1, then the third term in Equation (25) 
is negligible compared to the previous two. Hence, if the assumption 
is made that the frequency is high along with small perturbations, 
the simple wave equation with variable sound speed results: 


1 3 2 P 

a^(Y) 7t? 


V 2 P = Q 


( 26 ) 


where the source term 0 is represented by 


6(r) 

Q = q £(z-s)e iut 

nr 


a point source at z = s and r = 0, and from Equation (3) 


2, v 2 , AT 1 

a (z) = a (1 + — ) . 

M ® T® I + aZ 


At the surface (2=0) the boundary condition is 

- Z w = P (28) 

where w is the vertical component of the acoustic fluid velocity ob- 
tained from 3w/3t = - l/p(3P/3z) and Z is the normal acoustical impe- 
dance of the surface. Also, as z-*® a radiation condition requiring all 
waves to be propagating upward must be applied. 

The solution is obtained by first separating out the time depend- 
ance by defining a function G (z,r) such that 


P (z,r,t) = e iu)t G (z,r) . 


Taking the Hankel transformation (two-dimensional Fourier transform) 
of the resulting governing equation and boundary condition then yields 


d 2 G (u 2 

~~J + ( ~T 

A-t*- 


1 + aZ o 

- 6 2 )G 

1 + aZ + A i 


5 (Z-S) 


WPq Z = 0 


on z = 0. The function G(z,s) is the Hankel transform of G(r,z) and 
is defined by 
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G(z,B) = j 

G(z,r) r J 0 (Br) dr . 

(32) 

0 



The inverse 

Hankel transform yields G(r,z) 

and is gi ven by 

o (z,r) = / 

00 

G(z.b) B J 0 (er) dB 

■\ 

(33) 


Nayfeh [16] gives a method of obtaining an approximate solution 
of Equation (30) which is valid when X = tu/ ( a^a) >> 1. This requirement 
is the same as that found to be valid for Equation (26), so no new ap- 
proximation is required. The approximate solutions can be expressed as 


A 8 

G = i — — Mn(z)) + n — h 2 (n(z)) 
J 9 (z) (z) 


where 


3 2/3 , % 

n(z) = (— X) g(z) 


. . . 2 1/3 1/2 (1) 2 3/2 

1 (5) ■ (y> 5 H 1/3 ( T 5 ) 


and 


2 1/3 1/2 (2) , 2 3/2 

h 2 ( 0 = ( T } « H i /3 ( T' > 


(34) 


(35) 

(36) 


(37) 


Values of ine modified one-third order Hankel functions, h] and h 2 , 
are tabulated and their properties are discussed in Reference [17]. 


The function 
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2 (1 - ( (a^/u p)2) (1 + az + AT/T) - (aT/T) 

3 “ v g(z) (1 + a z + AT/T) 


(38) 


giving 


g3/2(z, S ) =/(l - (-%) ? )(1 

V U) 


AT 2 AT AT 

+ aZ + — ) - — (1 + aZ + — ) 

T T T 


1 AT ] , , 1 + * x 

— — 7 In ( ) 

2 T \|1 - ((a./u> 8 / 1 - * 


(39) 


where 


♦ = 


a_ 2 


AT AT 


(1 - (~e) )(1 + az + -) - - 


a„ 2 AT 

\( (1 - (—8) )(1 + az + — ) 


(40) 


This solution is the same as for the plane wave case [11] provided that 
cos 8 in the plane wave case is replaced by (aa,/a>)8. Here 8 is the 
incidence angle of the plane wave at infinity and a is the speed of 
sound at infinity. This result, obtained from the Hankel transforma- 
tion, shows that the point source problem in fact can be interpreted 
as a summation of the plane wave solutions. Before we can proceed to 
the final solution, it is necessary to examine the nature of the plane 
wave propagation associated with the point source problem. 

In the plane wave problem all real incidence angles were possible, 
and cos 8 could vary between zero and one. In the point source problem 
the range is more limited. First consider the waves propagating down- 
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ward from the source located at height s. Wave that are directed down- 
ward at a sufficient angle will reach the ground as in the plane wave 
case. Plane wave rays directed downward at an infinite height with 
incidence angles between 90° and some value e 0 , or 0 < B < Bo = 
(u)/a») cos e 0 = (tD/a*) J l/O+AT/T) can pass through the source loca- 
tion and reach the ground where they are reflected (see Figure 2). The 
turning points of these waves are below the ground surface. The initial 
angle 0 O characterizes the waves whose turning point is at the ground 
level, z = 0. Waves in this range of initial angle and their reflec- 
tions from the ground constitute two groups of possible waves and one of 
the forms of the solution. 

A second group of waves which is initially directed downward at 
the source are rays with angles at infinity in the range between e 0 and 
e z or b 0 < 6 < B z = (u/aoo) cos e 2 = ( w/a «) ( 1 + ctz)/(l + az + aT/T)). 
The angle 6 Z characterizes the wave whose turning point is located at 
the observer height z. This range of angles or of betas then character- 
izes waves that pass through the source location and whose turning 
point is above the ground level but below the observer (see Figure 3). 
These waves, before pass ng through the turning point, and after pass- 
ing through the turning point, contribute a second form to the solution. 

A third group of waves propagating downward at the source has 
angles at infinity in the range between e z and 0 S or Bz < B < Bs = 
(w/a*) cos 0 5 = (a)/a 00 )7((l + as ) / ( 1 + as + AT/T)). The value 0 S char- 
acterizes the wave with its turning point at the height s. It is the 
wave that leaves the source horizontally and is the limiting case of 
waves leaving the source and propagating downward. Waves in this 



Ray d 
the p 
refle 






Figure 3. Ray diagram for a plane wave passing through the point source location 
and being refracted before reaching the ground but below the receiver. 
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group have their turning points above the observer (see Figure 4), and 
thus cannot contribute an oscillating term to the total solution but 
rather have an exponential behavior. This behavior is further discussed 
in connection to the Saddle Point method of evaluation. This group of 
waves and their continuation after being refracted upward constitute 
a third form of the solution. 

The fourth form of the solution is due to waves which leave the 
source and propagate upward. These waves as shown in Figures 5 and 6 
have 9 s < 9 « 90° or 0 <0 < 0 S but are not of the same form as described 
earlier as their amplitude must be corrected for waves leaving the source 
rather than waves reflected or refracted upward. These waves constitute 
a fourth form of the solution. 

All of the waves of the four types described above can be seen 
in the ray diagram of the point source Figure 1, but a number of addi- 
tional types of waves need to be included in the total solution which 
do not appear in the ray diagram. The fifth form of the solution has 
6 S > 6 > C° or b s < B < u/aoo. These are waves with their turning points 
above the source and do not appear in the ray diagram but contribute to 
the full solution. In this case g(s,e) is a complex number with phase of 
tt/ 3 and the solution has an exponential behavior. 

The actual forms of the solution are obtained by requiring Equation 
(30) to satisfy the surface boundary condition, (31), at z = 0: at 
z = s, the source height, the solution must satisfy the condition that 

lim [G(6,z)] = lim [G( S , z ) ] (41) 

z+s_ z-s+ 

which requires continuity of pressure and the condition that 



Figure 4. Ray diagram for a plane wave passing through the point 
souce location and being refracted before reaching the 
ground and above the receiver. 




Figure 5. Ray diagram for a reflected plane wave passing through 
the point source and propagated upward. 




Figure 6. Ray diagram for a refracted wave passing through 
the point source and propagated upward. 
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r 3 G 3G q 

lim — - 1 im — = — 

z>s_ 13 zJ z+s + l_3zJ 2* 


( 42 ) 


which Is obtained by Integrating Equation (30) from z - s - t to z - 
s + c and taking the limits of t -*•0. At the turning point height 
the condition 

lim [G ( 6 , z ) ] = lim [G(B,z)] (43) 

Z+ZJP+ z*ZTP- 


is required if the turning point is above ground surface, z = 0. 
With some mathematical manipulation [li], the resulting solution can 
be written in the following form: 

For z > s. 


G = K h 2 (r,(B,z)) [ h! (n(B,s)) + R 0 (n(B,0)) h 2 (n(B,s)) ] (44) 


for w/a*, > B z > B s > B 0 > B > 0, 


G = K ( h 2 (n(B,z)) [h] (n(B,s)) 

. 2tt 

+ R 0 (n(B,0)e 13 ) h 2 (n(B,s))] 


(45) 


for u/a* > B z > B s > B > B 0 > 0* 


2t 

iT 


G = K ^ h 2 (n(B,z)) [h 1 (n(B,s)e J ) 

?v 2 t 

+ R„ (n(6,n)e' T ) h 2 (n(s.s)e’ 3 )] 


(46) 


for ui/d a >B z > B>Bs > Bo > 0 and 

r -h ■— 

G = K | h 2 ( n (6,z)e 13 ) [hj (n(B,s)e 13 ) 


2tt 


2t i 


* R 0 (n(e.0)e'^~) (n(B.s)e’^ )] 


( 47 ) 
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for 0 )/dao > B > B z > B s > B 0 > 0. For z < s 


G = K h 2 (n(B,s)) !>i (n(B,z)) + R 0 (n( B.0) ) h 2 ( n( B.z) )3 \ (48) 


for w/a*. > 0 S > 0 Z > B 0 > B > 0, 


2it 


G = K h 2 (n(B,s)) [h 1 (n(B,z)) + R 0 (n(B,o)e 13 ) h 2 (n(B,z))] 


49) 


for tu/ac > B s > B z > B > B 0 > 0, 


2tt 


G = K ^ h 2 (n(B,s)) [hj (n(B.z)e J ) 

2tt 2tt 

+ R 0 (n(B.O)e 1 3 ) h 2 (n(B,z)e 13 )] 


(50) 


for cu/ a*, > B s > B > B z > B 0 > 0 and 

2tt 2tt 

G = K( h 2 (n(6,s)e 13 ) [hj (n(8,z)e 13 ) 


2tt 


2tt 


+ R 0 (n(e,o)e 15_ ) h 2 ( n (g,z)e 13 )] 


for w/a*, > B > B s > B z > B 0 > 0, where 


(51) 


K = 


1 2i x“ 7 ^g z (B,s) y/g z (B,z) 


(52) 


R 0 - - 


T h! (n( B.0) ) + 1 ip r>l' ( n( 0,0) ) 
t h 2 (n( B.0) ) + i h ?‘ ( n( B,0) ) 


(53) 


i Z g ( 6,0 ) 
X = a X - - 


2 p3oo g z (b.o) 


(54) 


Z 2/3 32/3 

* = — X g z (B,0)(-) 

pa c 


( 55 ) 
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Note that g z (z) here denotes the derivative is taken with respect of 
az rather than z as in g' (:). These eight forms account for the five 
physical situations described above since s^me forms are repeated. These 
forms of solutions are useful in the way that the first term always 
represents the direct wave and the second the reflected or refracted 
wave. Tnis then completes the formal solution to the problem with the 
solution given in the form of an integral. Equation (33), the inverse 
transformation of the function G(B,z) given tj Equation (44) through 
(51). This is a very complex integral and approximations are necessary 
to proceed toward its numerical evaluation. 


CHAPTER IV 


EVALUATION OF THE FORMAL SOLUTION 

As we have seen in Chapter III, the problem that needs to be solved 
is the inverse Hankel transformation of the functions defined by 

G (r,z) = f G ( 8, z ) 6 J 0 ( Br) de (56) 

J o 

where G ( 8 , z ) is given in Equations (- 14 ) through (51) depending on the 
location of the receiver. Because of the complexity of the functions 
involved, it is unlikely that an exact solution can be obtained through 
an analytical approach. The second alternative that came into consi- 
deration was the direct numerical integration. However, in doing so, 
it woulo require the integral in the inverse Hankel transform to be 
evaluated at each point of interest. Of course, this could be done 

except that it requires a tremendous amount of computing time and, 
also, the understanding of the physical nature of the solution is 
reduced. The most satisfying remedy to this problem is to employ an 
approximate integration. The Saddle Point method, or the method of 
steepest descent, has been used successfully in the isothermal problem 
and appears to be useful here. 

The Saddle point method is a method of approximately integrating 


expressions of the form 
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a 

J 


c# . Cf (B) 

F ( 6 )e dB 


(57) 


where £ is a large parameter. If we interpret this integral as an 
integral in complex B space, the path of integration can be modified 
according to the rules of contour integration in the complex plane. 
What would apoear to be desired, then, is to find a point (or points) 
in the complex B plan? where the integrand of Equation (57) is a maximum 
and to distort the oath of integration to pass through tnis point in 
such a way that the exponential term in Eauation (57) decreases rapidly 
as one moves along the modified path of integration away from the 
point of maximum integrand. In the complex plane, however, a regular 
function does not have a maximum: but rather a point where 


3 f 
3B 


= 0 


B 



(58) 


is a saddlp point of the function, a point where the function surface 
looks 1 i k e a mountain pass, col, or saddle. This point is found to give 
the major contribution to the integral when the path is modified to give 
a rapidly decreasing integrand as one moves away from it. Approximately 
integrating Equation (57) then can be shown to give 


2tt 


(~ £f"(Bsp) 


f ( Bsp )e 5f(8 -» 


+ 1 O 


(59) 


wnere o is the angle at which the path of integration crosses through 
the saddle point. This method is described in more detail and with 
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more rigor in most texts on advanced mathematical methods. Reference 
[13], for example. This method is not a fool-proof one however, as 
the modified path of integration may intersect branch lines, pass around 
poles, or lead into other difficulties. Keller [18], however, asserts 
that the major part of the solution is obtained. Also if no saddle 
point occurs then the integral can be approximated as zero. Using this 
rather crude approach the integrals can be approximated relatively 
easi ly . 

In order to ultilize the Saddle Point method, some modifications 
of the integral are necessary. The method requires an integral involv- 
ing an exponential term to be evaluated from negative infinity to posi- 
tive infinity. This can be achieved by first noting that the Bessel 
function J 0 (pr) in Equation (56) can be written [19] as 

1 ( 1 ) ( 2 ) 

J 0 (BO = ~ [ H V '(er) + n V) ] ( 60 ) 
2 o o 


also 


H (1) (- Sr) = - H (2) (er) (61) 

o o 

and G(e,z) is an even function of e. Thus Equation (56) can be written 

as 



( 2 ) 

G (B.z) H (pr) de 
o 


(62) 


The exponential forms are obtained by then assuming the arguments of 
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all the Hankel functions involved are sufficiently large to allow re- 
placement by their asymptotic forms [19]: carrying out this asymptotic 
expansion process the resulting expressions are of the form (59) with 
the integrand of each integral consisting of two terms: one will be 
seen to represent the direct waves, and the second the reflected and/or 
refracted waves. This becomes clear when the derivative with respect to 
B is taken of the argument of the exponential terms to find the saddle 
points. It is found that 


3/2 

3(A g (B,z)) 
98 


= - F( 0 ,z) 


(63) 


where F(B,z) is the same function that occures in the ray tracing pro- 
ded ihit cos 9 is replaced by (ao«>/aj)6 where e is the incidence angle 
of the plane wave at infinity. Thus, the finding of the saddle points, 
which could be numerically quite difficult, has a clear physical inter- 
pretation in terms of the rays of the point source. For an observer 
located at a point in physical, (r,z), space which is "illuminated" 
by the point source, the saddle points of the approximate integrals 
found for the solution correspond to the B values of the two rays 
passing through that point as shown in Figure 1. Physically this is 
interpreted as indicating that the wave front travelling along the ray 
that passes through the source and receiver location is the major 
contributor to the sound levol that is experienced. 

At any point in the "illuminated" region of physical space the 
solution is then approximated, as described above, by two of the terms 
listed below. For direct waves, with z > s, region 1 and 2 of Figure 1, 
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g d = 


\f2^ K 2 e-MB^r + t } ) 
♦l' , ( 8 S p) 


(64) 


where 


*,(8) = -X g2/3 ( 6 , s ) + \ g 3/2 (s>z) 

2 


(65) 


and 0 sp is the root of 


r + (e sp ) = o . 


( 66 ) 


For direct waves with 0 < z < s, regions 3 sad 4 of Figure 1 


G d = 2rr K p e~ l(S sp * ^2 > 

f- ♦” ( 6 5 o ) 

where 

<P 2 (B) = X g 3/2 (b,s)- x g 3 / 2 (e,z)-l 

2 

and 0 S p is the root of 
r + $2 (B S p) = 0 • 

For reflected waves, regions 1,3, and 5 of Figure 1 


(67) 


( 68 ) 


( 69 ) 
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Gr = 


2* K 2 R 2 (e S p) e _i(0 sp r+, ^3) 

^-"TJTBspT 


where 


♦ 3 ( 6 ) = a g j/2 ( 6 * s ) + X g 3/2 (B,z) - 2x g 3/2 (B sp ,0) - ^ 


and B sp is the root of 
r + <f>2 ( 0 sp) = 0 • 


For refractPd waves which have not contacted thp caustic. Regions 
and 6 of Figure 1 


Gr = 


fa K 2 e- i(6 sp rt »3) 

/«>" (B ) 
v 4 *P 


R l< B sp>) 


where 


1 1 TT 


^4(6) = X g 3/? ( B , s ) + x g 3/2 (b.z) + — 


and 8 sd is the root of 


r + ^ (B sp ) = 0 . 


For refracted waves that have contacted the caustic, regions 


(70) 

(71) 

(72) 
4,5, 

(73) 

(74) 

(75) 
2,4. 


and 6 of Figure 1 
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, 2,V 2 e-'fV^s) 

g r = 

where 

it 

<t>5 = <t»4 “ ~ 


( 76 ) 


(77) 


and e S p is also a root (the smaller or the two roots for region 
6) of Equation (75). 

The derivative of g 3/2 ( 2 ,g) is defined by Equation (39). 


k 2 (*sp> 


_q__ 1 1- 1 

12tri \l it \Jg z ( e S p *z)^9 Z ( Bsp » s ) g 1/4 (Bsp< z ) 


] 

g 1/4 (B$ p ,s) 


(78) 


and 



C9) 



i + ip i (j\) 

T " ^ 1 (|x) 


1/3 

173 


V g(e,o) 

\ r g(M) 


(80) 


where R 0 , t and i|> are defined in Equations (53), (54) and (55). Thus, 

the problem has resolved itself into a relatively simple expression, 

but two difficulties remain. First, the expressions given above are not 
continuous as the receiver is moved from region to region in physical 
space. Secondly, no expression has been obtained that is valid in 

the shadow, region 7 in Figure 1. The first difficulty is easily 
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resolved by an ad hoc approach. The discontinuity results from the 
fact that on the boundary between regions the arguments of one of the 
modified Hankel functions is zero and it has been assumed that all 
of the arguments are largf enough to allow use of the asymptotic 
forms. The actual integrals are continuous at these boundaries but 
the asymptotic forms are not. Thus, by replacing the asymptotic forms 
of the modified Hankel functions in Equations (64), (67), (70), (73) 
and (76) by the modified Hankel functions themselves one obtains: in 
place of Equation (64) 


_ n[2tt" K 3 -i (B sp r - tt/2) 

Gd = 7~ o ~ — e h 2 (n(Bsp» z )) h l (n(B sp ,s)) 




iL , , 

3B* 1 


(81) 


for direct wave at z < s: in place of Equation (67) 


'[ 2 ir K 3 -i (B sp r - tt/ 2 ) 

Gq = e h j (n(Bsp.z)) h 2 (n(B S p,s)) 

32 / \ 

W '♦*’ 


(82) 


for direct waves at z > s: in place of Equation (70) 


Gr = 


>f?7 K 3 R 0 (B$p) -i ( Bsp r - n/2) 




3B 




hp ( n( Bsp,z) ) h 2 (n(Bsp» s )) 

(83) 


for reflected wave: in place of Equation (73) 


'J 2ir K3 Rj (B$p) -i ( Bsp r - u) 

Gr = f— : e h 2 (n(B S p* z )) h 2 (n(B S p.s)) 


98 


(^ 4 ) 


(84) 


for refracted waves that have contacted the caustic; and, in place of 
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Equation (76) 

yfz* Rj (8$p) r ” 

Gr = - - - e h 2 ( n( B S p * z ) ) h 2 ( &sp » s ) ) 

for refracted wave that have contacted the caustic. Here 


<3 


q 1 1 

24 i X 2/3 / 9 Z ( B sp ,z) /g z (B sp ,s) 


2B 

S P . 

irr 


( 86 ) 


The second problem, an expression valid in region 7, the shadow, 
can also be circumvented. Examination of Equations (84) and (85) 
the expressions valid in region 6 shows that these expression become 
in finite at the caustic, where 


3 2 4>4 _ 3 2 <{>5 

3B 2 38 2 


0 . 


(87) 


This result is easily understood by realizing that the refracted rays 
are given by r = F(s,b) + F ( z , 6 ) , the same as Equation (75) expressed 
in F. Recall that the caustic is defined by the largest radius that 
can be obtained for fixed observer height, z, and source height, s, or 
where 3r/3B = 0, thus yielding the above expression. This singularity, 
in the mathematical sense, results from expanding the argument of the 
exponential terms in the saddle point method in a Taylor series but 
retaining only second terms, which are zero in this case. In the 
physical sense it arises from the ray tube area going to zero and con- 
servation of energy. A result valid near the caustic has been obtained 
by Sachs and Silbiger [20]. Using their approach, which essentially 
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involves expanding the exponential terms about the ray tangent to the 
caustic rather than the saddle point, yields the expression 

G R ■ 2 M 6 c> <-- — > 1/2 e' 1( c r • T> Ai (( — - — ) 1/3 Ar) . 

( 88 ) 

h 2 (n(Bc» 2 )) h 2 ( n( 6 C • s ) ) 

where B c is the root of 

(s c ) = 0 (89) 

4 

which is valid near the caustic and on both the "illuminated" and the 
shadow side of it. Here 

Ar = r - r c (90) 

where the caustic radius 

r c = - 4> ' ( Be ) • (91 ) 

Tne applicability of Equation (88) is restricted not only to large 
frequencies, as we pointed out in the beginning of the discussion, 
but also to the vicinity of the caustic. Above the caustic, the 
range where Equation (88) is applicable is determined by a numerical 
matching between Equation (88) and the governing equations for the 
adjacent regions. In the shadow region, the criteria of the range of 
validity are not quite clear at the present time and will be left as a 
subject for the future study. 



CHAPTER V 


PROGRAM CONTENT 

With the results developed in the previous sections, we can now 
proceed to implement them with a computer program. It is hoped that 
this will enable us to examine thoroughly the behavior of sound inten- 
sity under the effect of thermal inhomogeneity, and as a result provide 
a better prediction scheme in the future. This program carries out 
the necessary procedures in the appropriate order, and its contents 
are described in the following. 

The first step is to determine at which region the acoustic re- 
ceiver resides. This is done by first evaluating the corresponding 
boundaries at the height of the receiver. There are four boundaries 
that separate the physical space into seven different regions, with 
each region characteri zed by two different rays as shown in Figure 1. 
Their governing equations are given as follows: 

Boundary I: Bj = F(z,0) . (92) 

This is the ray whose initial angle is zero, or at an angle parallel 
to the horizontal . 

Boundary II: Bjj = F(z,0s^) + F(s,0s^) . (93) 

This is the ray with its turning point at the ground (z = 0). Here 
0c corresponds to the initial angle at the source. 0c is obtained 
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from the fact that at the turning point 
Az + B = .0 


(94) 


or 


B = 0 


(95) 


where A and B are defined as before in Equation (10). Carrying out 
the mathematics we find that the angle 


Cos 9$ = 


1 + as + AT/T 


11 \j(l + AT/T)(1 + as) 


(96) 


or in terms of B 


0 ) 



1 

I + AT/T 


(97) 


Boundary lilt Bjjt = F ( s , 0$ j j j ( 1 2 ) ) * (98) 

This boundary does not represent one single sound ray but rather tne 
locus of all the turning points between the source and the ground. 
This is determined by the same equation that governs boundary II 


Az + B = 0 


(94) 


except that z now is varied depending upon the height in interest. 
Carrying out the appropriate steps gives 


Cos 0c (z) 
5 III 



1 + aZ 


+ aZ + AT/T 


1 + as + AT/T 


1 + as 


( 99 ) 


the wave, having a turning point above the receiver which we have ig- 
nored. Nevertheless, because of the reflective index R 0 (b), which in- 
involves complex value g^/^(z,B), its derivatives have to be eval- 
uated properly. Besides, it is very likely that this will be found to 
be useful if any further study of the approximation scheme is desired. 
As a result of choosing this particular branch cut, /-I is now being 
taken as -i and ln(-l) is equal to -1* instead of in. In terms of Car- 
tesian coordinates, the branch cut has shifted from the negative X-axis 
to the positive Y-axis for this purpose. All these modifications have 
to be taken into consideration in the programming process. 

Table 1 summarizes all the subroutine and functions being facil- 
itated in this program. 
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TABLE 1 


Subroutines and Functions 



SUBROUTINES 

Name 

Function 

FINREG 

Finding the region at which the receiver locates. 

FINDIR 

Finding the saddle point that corresponds to 
the direct wave. 

FINRFL 

Finding the saddle point that correspond to 
the reflected wave. 

FINRFC 

Finding the saddle point that corresponds to 
the refracted wave. 

INTEG 

Evaluating the pressure at the receiver's lo- 
cation. 

EVALB4 

Evaluating the location of boundary. 

ROOT 

Root finding routine ij^ng the secant method. 

HKFCT 

Evaluating tin modified Hankel functions. 


FUNCTIONS 


Name 


Function 



FNRAY 

632 

G32PI 

G32PI2 

632 PI 3 

GPIZ 

GPIZ2 

CG32 

CGPIZ 


Evaluating tne F^z.9 s ) function, 
evaluating the g3/2(z,B) function. 

Evaluating the first derivative of g3/2(z,g) with 


respect to 8. 

Evaluating :• 1 second derivative of g3/2(z,g) 
with respe^ "j S. 

Evaluating the third derivative of g3/2(z,e) 
with respect to 6. 

Evaluating the first derivative of 9^/2 (z,b) 
with respect to a. 

Evaluating the second derivative of g3/2(z,g) 
with respect to az. 

Complex version of G32. 

Complex version of GPIZ. 


CHAPTER VI 


DISCUSSION OF RESULTS 

Numerical evaluation «f the solution for the propagation of sound 
waves in a thermally stratified medium has been performed on a digitial 
computer. Plots of results can be found at the end of this section. 
Surface impedance values given by Chessel 's relationship [5] for grass 
have been used, with specific flow resistance chosen as 300 cgs units. 
Speed of sound high above the ground (z— ®) is chosen as 340 m/s, the 
standard atmospheric condition. The sound source has been chosen to 
be 3 meters above the surface with an arbitrary strength of unity for 
simplicity. Thesp values remain constant throughout the entire cal- 
culation. The various parameters considered for examination of their 
effects on sound propagation are tne angular frequency w, the tempera- 
ture parameters a and AT/T, which correspond to the vertical scale of 
temperature, respectively, m is chosen to be in the range between 
10,000 and 15,000 rad/sec. Using data obtained by NASA [21] and by 
Butterworth [22] and fitting these data to the temperature profile 
described in Equation (1), typical values of a are found to range 
from 1.2 to 5.9 nr 1 and aT/T from 0.0032 to 0.342 for the lapse condi- 
tions. The values used for examinations c Q based on the above range 
of variation. 

Figure 1 depicts the formation of regions in physical space with 
a = 2.5 m" 1 , AT/T = 0.025, and S = 3 m. With this figure as the refer- 
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ent, we ran see from Figures 7 through 10 that the variation of the tem- 
perature parameters a and aT/T resulted in shifting the relative loca- 
tions of the regions. This may result in switching from one governing 
equation to another as discussed in Chapter III but appears to 
have little effect on the general nature of the solution. 

Figure 10 through 20 are plots of sound pressure level (with an 
arbitrary reference) versus elevation at radii of 20, 40 and 80 m. 
Figures 10 to 12 are intended as the referent as parameters are varied 
in other drawings: they are for values of aT/T = 0.25, a =2.5 m _1 , S = 
3 m, the same conditions as the ray diagram of Figure 1, and frequency 
of about 1600 Hz (w = 10,000 rad/sec.) 

At r 20 m. Figure 10 shows a r ather normal looking interference 
Dattern with a mean sound level of about 32 dB. This radius is suffi- 
ciently close to the source that most of the acoustic energy is con- 
tributed from the direct and reflected waves: only waves leaving the 
source at very shallow angles are refracted upward before reaching the 
ground (see Figure 1). The regions of Figure 1 are indicated on Fig- 
ure 10 for reference. The oscillations, in this interference pattern, 
slowly increase in amplitude as it moves away f rom the source height. 
This is oue to the behavior of the strength of the reflected wave as 
shown in Figure 21. One conspicuous feature noticed in Figure 10 is 
the dramatic drop of sound level occurs at location near the source 
height S = 3 m: this violates our intuition that a higher sound level 
at that location should be experienced instead. This peculiar behavior 
near the source height and at small radii turned out to be the charac- 
teristic of the solution. Carefully examining the ray diagram in Figure 
22, which is plotted with the rays emitted in equal angular interval. 



Figure 7. Ray tracing diagram for 
and aP/T = 0.0?5. 
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Figure 9. Ray tracing diagram for 
and aT/T = 0.025. 



SOUND PRESSURE LEVEL VS HEIGHT 



SOUND PRESSURE LEVEL (dB) 





SOUND PRESSURE LEVEL VS HEIGHT 





SOUND PRESSURE LEVEL VS HEIGHT 



AT/T = 0.025, and w = 10,000 rad/sec 








SOUND PRESSURE LEVEL VS HEIGHT 



AT/T = 0.0?5, and oj = 15,000 rad/sec 






SOUND PRESSURE LEVEL VS HEIGHT 



Figure 15. SPL vs HEIGHT at r = 40 m, a = 2.5 
AT/T = 0.025, and oj = 12,000 rad/sec 







PRESSURE LEVEL VS HEIGHT 



rad/sec 










SOUND PRESSURE LEVEL VS HEIGHT 












B.O 








♦ 


66 



we noticed that the perpendicular distance between rays is further 
apart at this location than anywhere else at the same radius. If the 
equal acoustic energy is assumed to be present between rays, a location 
with a large distance between rays should be a region of low acoustic 
pressure. With this argument, the low sound level which occurs in 
the locations near the source height can then be accounted for. 

Figure 11 is a similar plot at r = 40 m. At this radius, the 
receiver has traversed a number of different regions and this accounts 
for the irregular interference pattern in the result. The mean sound 
level in this interference is about 26 dB, 6 dB less than the mean on 
the previous plot, as is expected. Near the source height a large 
drop of pressure takes place due to the ray spreading discussed above. 
One interesting feature observed from this plot is that the peak pres- 
sure occurs not on the caustic itself, as would be expected from the 
ray acoustics, but rather at some distance above it. In the shadow 
zone, the field is decaying exponentially with increasing distance 
from the caustic boundary. This is due to the diffraction effects 
which occur anytime a wavefront cross section is limited in some manner. 
In this case it is limited by an inability to penetrate the caustic 
boundary. This phenomenon is explained by the Huygen's principle 
discussed in some acoustic texts. Reference [15], for example. This 
particular location of the peak pressure and the decaying behavior 
can be seen in all the plots where there is a shadow region. At a 
height of about 1.6 m a discontinuity occurs in the sound for the 
integral level. This results from the higher approximation we used in 
the vicinity of the caustic boundary than the one for other locations. 

Figure 12 is plotted at r = 80 m. Note that the height scale 


V ' - 



V 
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does not continue to the ground at z = 0. It is because the approx- 


imation, Equation (88), we use in the vicinity of the causic is not 


valid for the entire shadow region as discussed in Chapter IV. Again 


the regions of Figure 1 are shown on the plot. Here, discontinuities 


in slope occur at about 4.2 and 2.2 m. The first is at the boundary 


between the direct and the refracted waves, resulting from the approx- 


imation error, and the second is accounted for by the numerical match- 


ing which has deliberately forced the solutions to be continuous 


in magnitude but not in slope. Again the mean sound level has de- 


creased about 6 dB to near 20 dB with the doubling of the distance 


from the source as compared to the previous figure. 


Figures 13 to 18 are plotted with the same parameters as the 


referent figures exceot that higher frequencies at u = 12,000 rad/sec 


(f = 1,910 Hz) and w = 15,000 rad/sec (f = 2,390 Hz) have been used. 


From these results, we found the increases in the frequency of the sound 


source have little effect on the general behavior of the sound level 


rather than an increased oscillation in the interference pattern. 


Plots with different temperature parameters a and AT/T can also be 


found in Figures 19 and 20. 


Figure 23 is a plot of sound level versus horizontal distance 


from the source. Again the regions of Figure 1 are shown along with 


the empirical model of Wiener and Keast [23]. The Wiener and Keast 


model suggests that a simple source model (corrected for atmospheric 


absorption which is not present here) yields the sound level between 


the source and shadow boundary. A simple source model represents a 


direct wave propagating out from the source in a homogenous medium 


with its attenuation accouting only for the geometrical spreading. 
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Thus, no interference pattern can be produced In this model but only 
the local mean sound level. Beyond the shadow boundary an empirically 
obtained excess attentuation is applied to the simple source model. 
This empirical model, based on the experimental data collected by Wiener 
and Keast, shows a good agreement with the local mean sound level in 
the interference region and good agreement with the slope of the 
decay of the s^und level in the shadow region. The Wiener and Keast 
model would be in better agreement with the present model in the shadow 
region if we applied the excess attenuation at the radial location 
where the sound level starts falling below the simple source value 
rather than at the actual shadow boundary where the sound level is 
always greater than the simple source. Figure 24 compares the excess 
attenuation seen in this case with the model of Wiener and Keast. 
Note that if Wiener and Keast did not correctly determine the shadow 
location the "constant" difference between their model and the present 
model would be accounted for. Thus, the best experimental data cur- 
rently available, which are summarized in the Wiener and Keast model, 
appear to be in good agreement with the present model. 

There are some mathematical difficulties in the present model even 
though it gives remarkable results. These results arise from the asymp- 
totic behavior of the Hankel functions which restrict their arguments 
from being too large or when overflow occurs during calculation. In 
the present case, the arguments of the Hankel functions involve the 
g 3 /2 (z,B,) ari( j x, or (u/aa), terms; thus variation cT any of these vari- 
ables z, 6, u> , a, and a can result in too big a value for the Hankel 
function subroutine to handle, which in turn will appear as an overflow 
error in the computer program. As far as this program is concerned, the 
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upper limit for the argument of the Hankel function is a magnitude 
of about 27. The range of u and a has been discussed in a previous 
paragraph; and since a In practice does not vary significantly, the 
question remains on finding the range of the variables z and B for which 
the Hankel function can be evaluated properly. In doing so, it is 
necessary to first examine the behavior of the g3/2(z,g) function 
with respect to its arguments z and B. It is found that this function 
increases with increasing z and decreasing b and vice versa. Since z 
and B represent the height of receiver and the initial angle of the 
sound ray (B s - (w/a*) /(I + as)/(l + as + aT/T) cos 9 s ) that reaches 
the receiver, respectively, their limits actually form a boundary in 
the pnysical space beyond that which the Hankel function cannot be 
evaluated. For example, with x set equal to 11.8 rad and a height of 
5 meters above the ground, the minimum horizontal distance r from the 
sourd source is around 10 meters. At radii below this the problem 
will occur. The limit of r ratner than b is used because the saddle 
point B is totally depended on the location of the receiver specified 
by z and r, which are the input data for the computer program. On the 
other hand, in tracing the height along a constant radius of 40 meters, 
we found the maximum value of z is about 10 meters. Beyond that, an 
error will occur. Fortuitously, the locations corresponding to these 
large arguments are usually of no significant interet: they are places 
either in the close vicinity of the sound source or very high above it. 



CHAPTER VII 


CONCLUSION 

A model of acoustic prorogation in a thermally innomogeneous medium 
has been developed by Van Moorhem and numerically implemented on a Vax 
computer at the University of Utah. The temperature profile considered 
is a lapse condition, with temperature decreases in neight, along with 
a finite impedance ground surface using Vessel's model. The results 
appear to be in good agreement with the empirical model of Wiener and 
Keast, which by far. has the most reliable experimental data available. 
Yet further improvement is necessary for the present model. Poles that 
may occur in the integration process which are believed to lead to a 
surface wave term have not been considered here. Also the behavior of 
sound propagation in the shadow region far away from the caustic is yet 
to be examined. Investigation of these two subjects has been undertaken 
by Yiping Ma [24] and their resolutions are expected to be available 


in the near future. 


APPENDIX 


PROGRAM ACOUSTICS AND ITS USAGE 

This FORTRAN program carries out the numerical calculations as dis- 
cussed in the preceding sections. It was written on the VAX/VMS 

system at the Mechanical Engineering Department ot the University 

of Utah. With the input of the receiver's location along with the 

appropriate parameters, the program will evaluate the corresponding 
sound pressure level. Except for the location of the receiver, all 
the other parameters are assigned with default values which are com- 
monly used. The command 'namelist' is used to provide a convenient 

way to specify changes in the parameters of concern. The default 
values are assigned as follows 


Mathematical 

Variable 

FORTRAN 

Variable 

Default 

Value 

a 

SPEED 

340 m/s 

at/t 

DTOT 

0.025 

u> 

OMEGA 

10,000 Hz 

a 

ALPHA 

2.5 m-1 


In order to execute tne program, simply type in 'RUN ACOUSTIC' 
followed by the changes of defaulted parameters if desired. 

Example: $RUN ACOUSTIC 

SVALUES 


SEND 
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The above commands will initiate the execution of the program using 
the default values since no changes are specified. 

Example 2 : $RUN ACOUSTIC 

$VALUES 
OMEGA = 15000. 

ALPHA = 2.5 
SEND 

The above commands will initiate the execution of the program substi- 
tuting OMEGA = 15,000.0 and ALPHA = 3.0 for their default values. The 
program will then ask for the location of the receiver and the number 
of points that need to be evaluated. Here 'ELEV' is the height and 
' H D 1ST ' is the horizontal distance of the receiver's location. 1 R LOW * 
is the lower limit of the height to be evaluated and 'NO' is the number 
of values calculated between 'ELEV' and ' R L 0W . ' 

Example 3. ENTER INPUT ELEV. HDIST : 

5.0, 20.0 

ENTER INPUT RL0W, NO: 

0.0, 100 

These input values will enable the program to evaluate the receiver's 
location from 5.0 meters to ground (0.0 meter) along a constant radius 
of 20.0 metprs with intervals of 5.0/100 or 0.02 meters. If only one 
point in space is being evaluated, ' R LOW ' will be the same as 'ELEV' 
and NO wi 1 1 be 1 . 

All calculations are performed in double precision arithmetic. 

After the program has carried out the entire calculation, a result 
of sound pressure level versus height will store in a data file named 
' F OR 026 . ' A program named 'MGRAPH' is available and recommended for 
a quick view of thp plot of the result if the user runs this program 
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on tne same system as it the program was written. A copy of program 
ACOUSTIC can be found in the remainder of this appendix. 


OOOOOOCJOOO ooooooooooooooo o ooooooo 
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PROGRAM ACOUSTIC 


*** THIS PROGRAM DETERMINES THE SOUND PRESSURE AT ANY 
*** GIVEN LOCATION IN A THERMALLY STRATIFIED MEDIUM 
*** WITH THE TEMPERATURE PROFILE DISCUSSED IN CHAPTER 
*** II. A NUMBER OF LOCATIONS CAN ALSO BE EVALUATED 
*** IN THE EVENT THAT A PLOT OF SOUND PRESSURE LEVEL 
*** VERSUS HEIGHT IS DESIRED. 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*1 6 PDIR, PRFL, PRFC1 ,PRFC2,PCAU,PD(300),PR(300) , 
1 RG1 (300),RG2(300),PCU(300) 

DIMENSION A(300),SADRFC(2) 

REAL*8 MAG (300) 

INTEGER REGION 

COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT /TWO/S 
NAMELIST /VALUES/SPEED, OMEGA, ALPHA, DTOT, S 


*** VARIABLE USAGES: 


w w u 

Inrf 

*** 

•iHHt 

**•* 

*** 

*■#-* 


S = HEIGHT OF THE SOUND SOURCE'S LOCATION 
SPEED = SPEED OF SOUND OUTSIDE THE THERMALLY 
STRATIFIED REGION (AT INFINITE HEIGHT) 
OMEGA = ANGULAR FREQUENCY 
ALPHA = A TEMPERATURE PROFILE PARAMETER 
DTOT = (DELTA-T ) / (T-INFINITY ) 

ELEV = HEIGHT OF THE RECEIVER'S LOCATION 
HD13T = HORIZONTAL DISTANCE BETWEEN RECEIVER 
AND THE FOINT SOURCE 


S=3.0D+00 
SPEED=3*4D+02 
0ME)GA=1 .0D+04 
ALPHAS. 5D+00 
DT0T=2.5D-O2 

PI =4 • 0D+O0*DAT AN ( 1 . OD+OO ) 

READ (5, VALUES) 

PRINT*, 'ENTER INPUT ELEV,HDIST: ' 

READ*, ELEV, HDIST 

*** RLOW IS THE LOWER LIMIT OF THE HEIGHT TO BE *** 
*** EVALUATED. (FOR THE CASE THAT A NUMBER OF *** 
*** LOCATIONS NEED TO BE CALCULATED). *** 

*** NO IS THE NUMBER OF VALUES CALCULATED *** 

*** BETWEEN ELEV AND RLOW. *** 


*** 
M M W 

inrT 

*-*•* 
m m y 

*** 

*** 


*-*-* 

*** 
M M M 

ITU I 

*-*•* 

*-** 

*-*■* 


PRINT*, 'ENTER INPUT RLOW, NO: ’ 
READ*, RLOW, NO 


oooooooooooo o ooo o o oooo ooo 
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PRINT 1 

1 FORMAT (T) 

IF (NO .EQ. 1 )THEN 
RINC=0.0 
GO TO 3 
HJDIF 

*** RINC IS THE INCREMENT *** 

RINC= (ELEV-RLOW ) / (NO-1 ) 

ELEV=ELEV-f-RINC 

*** ICOMP IS A CONTROL STATEMENT FOR THE *** 
*** NUMERICAL MATCHING. *** 

IC0MP=C 


DO 200, K=1 ,N0 . 

ELEV=ELEV-RINC 
PRINT 5, ELEV,HDIST 

5 F0RMAT( 'O' , 'ELEV = ' ,F7-3,2X, ’HDIST = \F7-3) 

IF (ELEV .LT. 1.E-10)ELEV=0. 000001 

*** FUNDING THE REGION AT WHICH THE RECIVER LOCATES *** 
CALL FINRBG(S, ELEV, HDIST, REGION) 

GO TO (10, 20, 30, 55, 50, 55, 55), REGION 


*** THE SUBROUTINES ARE NAMED AS FOLLOWS: *** 

!l M M WWW 

*** FINDIR - FINDING THE SADDLE POINT FOR THE *** 

*** DIRECT WAVE *** 

*** FINRFL - FINDING THE SADDLE POINT FOR THE *** 

*** REFLECTED WAVE *** 

*** FDIRFC - FINDING THE SADDLE POINT FOR THE *** 

*** REFRACTED WAVE *** 


*** INTEG - EVALUATING THE SOUND PRESSURE WITH *** 
*** THE EQUATIONS DEVELOPED IN CH. Ill *** 


10 CALL FINDIR (S, ELEV, HDIST, SADDIR) 

IDBC=1 

CALL INTEG ( S , ELEV , HDIST , REGION, SADDIR , PDIR , IDEC ) 
PRINT* 

CALL FINRFL (S, ELEV, HDIST, SADRFL) 

IDBC=2 

CALL INTEG (S, ELEV, HDIST, REGION, SADRFL, PRFL, IDEC) 
GO TO 100 

20 CALL FINDIR (S, ELEV, HDIST, SADDIR) 

IDBC=1 

CALL INTEG ( S , ELEV , HDIST , REGION , SADDIR, PDIR, IDEC ) 
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PRINT* 

CALL FDJRFC(S, ELEV, HDIST, REGION, SADRFC) 

IDEC=2 

CALL INTEG(S,ELEV,HDIST, REGION, SADRFC(1 ),PRPC1 ,IDBC) 

C 

C *** NUMERICAL MATCHING FOR THE APPROXIMATED SOLUTIONS *** 

C *** BETWEEN REGION 2 AND THE SHADOW REGION *** 

C 

IF (REGION .EQ. 4) THEN 

ORIG=DSQRT (DREAL( PDIR+PRFC1 )**2+DIMAG(PDIR+PRFC1 )**2) 
RMOD=DSQRT ( DREAL ( PCAU ) **2+DIMAG ( PCAU ) **2 ) 
m PCT=ABS ( ORIG-RMOD ) /ORIG 
IF (ERPCT .LT. 0.02) THEN 
IDEC=3 
IC0MP=1 
ENDIF 
ENDIF 
GO TO 100 
C 

30 GO TO 10 
C 

40 GO TO 55 
C 

50 CALL FINRFL(S,ELEV,HDIST,SADRFL) 

IDEC=1 

CALL INTEG(S, ELEV,HDIST, REGION, SADRFL,PRFL,IDEC) 

PRINT* 

CALL FINRFC (S,ELEV,HDIST, REGION, SADRFC) 

IDEC=2 

CALL INTEG(S,ELEV,HDIST, REGION, SADRFC (1 ),PRFC1 ,IDEC) 

GO TO 100 
C 

C *** EVALB4 - EVALUATING BOUNDARY IV *** 

C 

55 CALL EVALB4(S,ELEV,ANGLE4,B0UND4) 

SADCAU= ( OMEGA/SPEED ) *DSQRT ( ( 1 . +ALPHA*S ) / ( 1 . +ALPHA*S+DT0T ) ) * 
1 DC0S(ANGLE4) 

PRTA=G32PI3 ( S , SADCAU ) 

PRTB=G32PI 3 ( ELEV , SADCAU ) 

ZETA=(2./(PRTA+FRTB) )**(1 ./3.) 

DELR=HDIST-B0UND4 

ARGU= ABS ( ZETA*DETR ) ** ( 3 • /2 . ) 

GO TO 70 
C 

60 CALL FLNRFC (S, ELEV, HOIST, REGION, SADRFC) 

IDEC=1 

CALL INTEG(S, ELEV, HDIST, REGION, SADRFC ( 2 ),PRFC1 ,IDEC) 

IDEC=2 

CALL INTEG(S, ELEV, HDIST, REGION, SADRFC(1 ),PRFC2,IDEC) 

C 

C *** NUMERICAL MATCHING FOR THE APPROXIMATED SOLUTIONS *** 

C *** BETWEEN REGION 6 AND THE SHADOW REGION *** 

C 


0RIG=DSQRT(DREAL(PRFC1+FRPC2)**2+DIMAG(FRFC1+PRFC2) # *2) 
RMOD=DSQRT ( DREAL ( PCAU ) **2+DIMAG ( PC AU ) **2 ) 

IRPCT=ABS ( ORIG-RMOD ) /ORIG 
IP (ERPCT .Iff. 0.02) THIN 
IDEC =3 
IC0MP=1 
INDIP 
GO TO 100 

IF ((REGION .EQ. 7) -AND. (ARGU .GT. 1.0)) THEN 
GO TO 350 
FTT^TF! 

IDEC=3 

CALL IN7EG(S,ELEV,HDIST, REGION, SADCAU,PCAU, IDEC) 

*FiN i yrp 

IP (ICOMP .EQ. 1)G0 TO 100 
IF (REGION .EQ. 4)G0 TO 20 
IP (REGION .EQ. 6)G0 TO 60 


A(K)=ELEV 
CT= 10000. 

PD(K)=PDIR*CT 
PR(K)=PRFL*CT 
RG1 (K)=PRPC1*CT 
RG2 (K ) =PRFC2*CT 
PCU(K)=PCAU*CT 

*** WRITING THE VALUES OF THE HEIGHT AND ITS *** 

*** CORRESPONDING SOUND PRESSURE LEVEL INTO *** 

*** the DATA PILE. *** 

GO TO (11 ,21 ,11 ,21 ,51 ,61 ,71 ), REGION 
MAG(K)=DSQRT(DREAL(PD(K)+PR(K) )**2+DIMAG(PD(K)+PR(K) ) 

** 2 ) 

GO TO 81 

IF (IDEC .EQ. 3)G0 TO 71 

MAG(K)=DSQRT(DREAL(PD(K)+RG1 (K))**2+DIMAG(PD(K)+RG1 (K)) 
** 2 ) 

GO TO 81 

MAG(K)=DSQRT(DREAL(PR(K)+RG1 (K))**2+DIMAG(PR(K)+RG1 (K)) 
** 2 ) 

GO TO 81 

IP (IDEC .EQ. 3)G0 TO 71 

MAG(K)=DSQRT(DREAL(RG1 (K)+RG2(K) )**2+DIMAG(RG1 (K)+RG2(K)) 

** 2 ) 

GO TO 81 

MAG ( K ) =DSQRT ( DREAL ( PCU ( K ) ) **2+DIMAG ( PCU ( K ) ) **2 ) 

WRITE (26,^1) 20*DLOG10(MAG(K)),A(K) 
P0RMAT(1X,F10.5,2X,P7.3) 


CONTINUE 
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350 STOP 
END 


SUBROUTINE FINREG ( S , ELEV , HDIST , RBCilQN ) 

*** FINDING THE REGION AT WHICH THE RECEIVER LOCATES *** 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

INTEGER REGION 

COMMON /ONE/SPEED , OMEGA , ALPHA , DTOT 

THE FOLLOWING STATEMENT IS USED BECAUSE OF THE 
PROBLEM ARISED FROM THE DOUBLE PRECISION. 

IF (ABS(ELEV-S) .LT. 0.1E-10)G0 TO 5 

IF (ELEV .LE. S) THEN 


IF (ELEV .LT. IE-5) GO TO 8 

*** ANGLE2 - THE INITIAL ANGLE OF THE RAY WHICH HAS *** 
*** IKE TURNING POINT AT GROUND (Z=0) 

*** B0UND2 - BOUNDARY II 

ANGLE2=DAC0S (DSQRT ( ( 1 . +ALPHA*S+DTOT ) / ( ( 1 . +DT0T ) 

1 *(1 .+ALPHA*S)))) 

B0UND2=FNRAY (ELEV , ANGLE2 ) +FNRAY ( S , ANGLE2 ) 

IF ((ELEV .GT. 2.99999) -AND. (ELEV .LT. 3*00001 ) )THEN 
B0UND3=0*0 
GO TO 10 
ENDIF 


W M W 

**-# 

*** 
n M w 

WWW 


ANGLE3 - THE INITAL ANGLE OF THE RAY WHICH HAS *** 

THE TURNING POINT AT THE RECEIVER'S *** 

HEIGHT *** 

B0UND3 - BOUNDARY III *** 


ANGLE3=DACCS ( DSQRT ( ( 1 . +ALTHA*S+DTOT ) * ( 1 . +ALPHA* 

1 ELEV)/((1 .+ALPHA*ELEV+DT0T)*(1 .+ALPHA*S)))) 

B0UND3=FNRAY ( S , ANGLE3 ) 

IF (ELEV .LT. 1E-5)B0UND2=B0UND3 
C 

10 IF ((HDIST .GT. B0UND2) .AND. (HDIST .GT. B0UND3)) THEN 

IF (ELEV .GT. 1E-5)THEN 

CALL EVALB4(S.ELEV,ANGLE4,B0UND4) 

ELSE 

B0UND4=B0UND2 

ENDIF 


ooo o ooooo ooo 
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I? (HDIST .LE. B0UND4) THEN 
REGI0N=6 

REGI0N=7 

ENDIF 

ELSE 

IF ((HDIST .GT. B0UND2) .AND. (HDIST .LT. B0UND3) ) THEN 
REGI0N=4 

EI£EIF ((HDIST .GT. BOUND?) .AND. (HDIST .LT. 

1 B0UND2) ) THEN 

REGI0N=5 
ELSE 
REGION =3 
ENDIP 
ENDIF 
C 
C 

ELSE 

C 

C 

ZERO=0.0 

*** B0UND1 - BOUNDARY I *** 

B0UND1 =FNRAY(ELEV, ZERO ) 

*** ANGLE2 - THE INITIAL ANGLE OP THE RAY WHICH HAS THE 
*** TURNING POINT AT GROUND (Z=0) *** 

*** B0UND2 - BOUNDARY II *** 

ANGLE2=DACOS (DSQRT ( ( 1 . +A1PHA*S+DT0T )/((!. +DTOT ) 

1 *(1 .+ALPHA*S)))) 

B0UND2=FNRAY ( ELEV , ANGLE2 ) +FNRAY ( S , ANGLE2 ) 

IF ((HDIST .LT. B0UND2) .AND. (HDIST .LT. B0UND1 ) ) THEN 
REGI0N=1 

EL3EIF ((HDIST .LT. BOUND 1 ) .AND. (HDIST .GT. B0UND2)) THEN 
REGI0N=2 

EISEIF ((HDIST .LT. B0UND2) .AND. (HDIST .GT. B0UND1 )) THEN 
REGIQN=5 
ELSE 

*** EVALB4 - EVALUATING THE BOUNDARY IV (THE CAUSTIC) *** 

CALL EVALB4(S,ELEV,ANGLE4,B0UND4) 

IF (HDIST .LT. B0UND4) THEN 
R 1 XjI0N=6 
GO TO 20 
ENDIF 
REGION =7 
ENDIF 
C 
C 


ooo ooooo 


82 


HJDIP 

C 

20 PRINT 25, REGION 

25 FORMAT (IX,' IT LOCATES AT REGION ',12) 

C 

RETURN 

END 


SUBROUTINE FINDIR(S,ELEV,HDIST f SADDIR) 

*** FINDING THE SADDLE POINT FOR THE DIRECT WAVE *** 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /ONE/SPEED , OMEGA , ALPHA , DTOT 
C 

THETA=DAT AN ( ABS ( ELEV-S ) /HDIST ) 

IF (ELEV .LT. S) THEN 

ANGLE3=DAC0S ( DSQRT ( ( 1 . +ALFHA*S+DTOT ) * ( 1 . +ALPHA*ELEV ) / 
1 ((l.+ALPHA*ELEV+DT0T)*(l.+ALPHA*3)))) 

ENDIF 

QUAN= ( OMEGA/SPEED ) *DSQRT ( ( 1 . +.ALPHA*S ) / 

1 ( 1 . +ALPHA*S+DTOT ) ) 

XL=QUAN*DCOS (THETA ) 

C 

IF (ELEV .GT. S) THEN 
XR=QUAN 

ELSEIF (THETA .GT. ANGLE3) THEN 
XR=XL-0.02 

P 

XL=QUAN*DCOS ( ANGLE3 ) 

XR--XL-0.001 

ENDIF 

C 

IC0UNT=0 

IF (ELEV .GT. S) THEN 
SI =-1.0 
ELSE 
SI=1 .0 
ENDIF 
C 

10 IC0UNT=IC0UNT+1 

IF (ICOUNT .GT. 100) THEN 

PRINT*, 'SOMETHING WRONG WITH FINDIR! ! ! ' 

GO TO 20 
ENDIF 

FXL=HDI ST+SI * ( G32PI ( S , XL ) -G32PI ( ELEV , XL ) ) 

FXR=HDI ST+SI * ( G32PI ( S , XR ) -G32PI ( ELEV , XR ) ) 

CALL ROOT(XL,XR,FXL,FXR, SADDIR) 

IF (SADDIR .EQ. 0.0)G0 TO 10 


PRINT 30,SADDIR 

FORMAT (IX,' SADDLE POINT FOR THE DIRECT WAVE IS ’,F7.4) 

PI =4 • OD+00*D ATAN ( 1 .OD+OO) 

OA =OMEGA/ SPEED 

ANGLE=DACOS ( SADDIR/ ( OA*DSQRT ( ( 1 . +ALPHA*S ) / ( 1 . +ALPHA*S+DTOT ) ) ) ) 
*180. /PI 
PRINT 40, ANGLE 

FORMAT (IX,' INITIAL ANGLE FROM THE SOURCE IS' F5.2,' DEG') 
CDIST=SI*(G32PI (ELEV, SADDIR)-G32PI (S, SADDIR) ) 

PRINT 50, CDIST 

FORMAT (IX,' CHECK: HDIST FOR THE DIRECT WAVE IS \F7.4) 

RETURN 

END 


SUBROUTINE FINRFL(S,ELEV,HDIST,SADRFL) 

*** FINDING THE SADDLE POI NT FOR THE REFLECTED WAVE ***■ 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /ONE/SPEED , OMEGA , ALPHA , DTOT 

PI =4 • 0!>+00*DATAN ( 1 . OB+OO ) 

ANGLE2=DAC0S (DSQRT ( ( 1 . 4ALPHA*S+DT0T ) / ( ( 1 . +DTOT ) * 
(1.+ALPHA*S)))) 

GUESS 1 =ANG hIP. 

GUESS2=Pl/2. 

I COUNT =0 

ZERO=0.0 

IC0UNT=IC0UNT+1 

IF (ICOUNT .GT. 500) THEN 

PRINT*, ’SOMETHING WRONG WITH FINRFL! ’ 

GO TO 40 
ENDIF 

REFL1 =HDI3T-FNRAY(ELEV, GUESS 1 )-FNRAY(S, GUESS 1 ) 

+2 . *FNRAY ( ZERO , GUE3S1 ) 

REFL2=HDIST-FNRAY ( ELEV , GUESS2 ) -FNRAY ( S , GUESS2 ) 

+2*FNRAY ( ZERO , GUES32 ) 

CALL ROOT ( GUESS 1 ,GUESS2,REFL1 ,REFL2, SADDLE) 

IF (SADDLE .EQ. 0.0)G0 TO 10 

SADRFL= ( OMEGA/SPEED ) *DSQRT' ( ( 1 . +ALPHA*S ) / ( 1 . +ALPHA*S+ 

DTOT) )*DCOS( SADDLE) 

PRINT 20.SADRFL 

FORMAT (IX,’ SADDLE POINT FOR THE REFLECTED WAVE IS ’,F7.4) 
PRINT 30,SADPLE*180./PI 

FORMAT (IX,’ INITIAL ANGLE FROM THE SOURCE IS \F5.2,’ DEG’) 


ooo ooooo 
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CDIST=FNRAY(ELEV , SADDLE )+FNRAY ( S , SADDLE )-2 . *FNRAY ( ZERO , SADDLE ) 
PRINT 35.CDIST 

35 F0RMAT(1X, 'CHECK: l^IST FOR THE REFLECTED WAVE IS ',F8.4) 

C 

40 RETURN 
END 


SUBROUTINE FINRFC (S, ELEV, HDIST, REGION ,SADRFC) 

*** FINDING THE SADDLE POINT FOR THE REFRACTED WAVE *** 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

INTEGER REGION 
DIMENSION SADRFC(2),RFC(2) 

COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT 

ANGLE2=DAC0S ( DSQRT ( ( 1 . +ALPHA*S+DTOT ) / ( ( 1 . +DTOT ) * 

1 (l.+ALPHA*S)))) 

IF (ELEV .LT. S) THEN 

ANGLE3=DAC0S ( DSQRT ( ( 1 . +ALPHA*S+DTOT ) * ( 1 . +ALPHA*ELEV ) / 
1 ( (1 .4-ALPIiA*ELEV+DT0T)*(1 .+ALPHA*S) ) ) ) 

ENDIF 
C 

NUM=1 

IF (REGION .EQ. 5) THEN 
GUESS1 =0.0000001 

IF (ELEV .LT. S) GUESS 1 =ANGLE3+. 0000001 
GUESS2=0.05 

10 REFL2=HDIST-FNRAY ( ELEV , GUESS2 ) -FNRAY ( S , GUESS2 ) 

IF (REFL2 .GT. 0.0) THE] 

GUESS2=GUESS2-fO . 02 
GO TO 10 
ENDIF 

ELSEIF (REGION .EQ. 4) THEN 
GUESS 1 =ANGLE2 
GUESS2=ANGLE3+. 000001 
EI^EIF (REGION .EQ. 2) THEN 
GUESS 1 =ANGLE2 
GUESS2=0. 0000001 
ELSEIF (REGION .EQ. 6) THEN 

CALL EVALB4(S,ELEV,ANGLS4,B0UND4) 

NUM=2 

ENDIF 

C 

DO 40,J=1 ,NUM 

C 

IF ((REGION .EQ. 6) .AND. (J .EQ. 1)) THEN 
GUES31 =0.0000001 

IF (ELEV .LT. S) GUESS 1 =ANGLE3+0. 00000001 


oooo ooooo 


85 


GUESS2=ANGLE4 

ELSEIP ( (REGION .EQ. 6) .AND. (J .EQ. NUM)) THEN 
GUESS1 =ANGLE4 
GUESS2=ANGLE2 
ENDIF 
C 

IC0UNT=0 

C 

50 ICCUNT=IC0UNT+1 

IP (ICOUNT .GT. 100) THEN 

PRINT*, 'SOMETHING WRONG WITH FINRFC ! ' 

GO TO 50 
ENDIF 
C 

FXL=HDIST-FNRAY(ELEV,GUESS1 ) -FNRAY(S, GUESS 1 ) 
FXR=HDIST-FNRAY ( ELEV , GUESS2 ) -FNRAY ( S , GUESS2 ) 

CALL ROOT ( GUESS 1 ,GUESS2,FXL,FXR,RFC(J) ) 

IF (RFC(J) .EQ. 0.0)G0 TO 30 
C 

SADRFC (J )= (OMEGA/SPEED )*DSQRT ( ( 1 .+ALPHA*S)/(1 .+ALPHA* 

1 S+DTOT ) ) *DCOS (RFC ( J ) ) 

C 

PRINT 35, SADRFC (J ) 

35 FORMAT (IX, 'SADDLE POINT FOR THE REFRACTED WAVE IS ',F7.4) 
ANGLE=RFC ( J )*1 80. /3 . 1 41 592654 
PRINT 37, ANGLE 

37 FORMAT (IX, 'INITIAL ANGLE FROM THE SOURCE IS ' ,F5.2, 'DEG' ) 
CHK=FNRAY(ELEV,RFC(J))+FNRAY(S,RFC(J ) ) 

PRINT 38,CHK 

38 FORMAT (IX, 'CHECK: HDIST FOR THE REFRACTED WAVE IS ',F8.4) 
C 

40 CONTINUE 
50 RETURN 
END 




SUBROUTINE INTEG ( S , ELEV , HDI ST , REGI ON , SADDLE , PINT , IDEC ) 

*** EVALUATING THE SOUND PRESSURE WITH EQUATIONS *** 
*** GIVEN IN SECTION 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*1 6 ZI.PRT1 ,K1 ,ARG1 ,ARG2,H1 ,H2,H1 P,H2P,PRT4A, 

1 PRT4B,PRT5A,PRT5B,PRT6,ARG3,T1 ,T2,R,PRT8, 

2 PRT6A , PRT6B , PINT , G32C , CMPSAD , CGPIZ , CGPIZ2 , 

3 ARG,AI,BI,AIP,BIP,PINT1 
REAL*8 LAND A, IM I 

LIT EGER REGION 

COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT 
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C 

f PI =4 . OD+CO*DATAN ( 1 . OD+OO ) 

SIGMA=300. 

*** ZI IS THE NORMAL IMPEDANCE FROM CHESSEL’S MODEL *** 

REI=1 .+9.08*(2.*PI*300./OMBGA)**0.75 
IMI=-11 .9*(2.*PI*300./0MBGA)**C.73 
ZI=DCMPLX(REI , IMI ) 

LANDA=OMEGA/(SPEED*ALTHA) 

ZER0=0.0 

CMPSAD=SADDLE*DCMPLX ( 1 . ,0.0) 

PRT1=DCMPLX( 0.0,1 .0)*24.*LANDA**(2./3. ) 

PRT2=DSQRT ( GPIZ ( ELEV , SADDLE ) *GPIZ ( S , SADDLE ) ) 

PRT3=DSQRT ( (2- ^SADDLE ) / ( PI *HDIST ) ) 

K1 =PRT3/(PRT1 *PRT2) 

ARG1 = ( 3 . *LANDA/2 . *G32 ( S , SADDLE ) ) ** ( 2 . /3 • ) 

ARG2= ( 3 • *LANDA/2 . «G32 ( ELEV , SADDLE ) ) ** ( 2 . /3 • ) 

CALL H(ARG1 ,H1 ,H2,H1P,H2P) 

PRT4A=H1 
PRT4B=H2 

CALL H(ARG2,H1 ,H2,H1P,H2P) 

PRT5A=H1 
PRT5B=H2 

PRT 6 =CDEXP ( DCMPLX ( 0 . 0 , -1 .0)*(SADDLE*HDIST-PI/2. ) ) 
FRT7A=G32PI2 ( S , SADDLE ) 

PRT7B=G32PI2(ELEV, SADDLE) 

IP ((IDEC .EQ. 2) .OR. (TDEC .EQ. 3)) GO TO 100 
GO TO (10,10,10,10,100,100) REGION 

10 IP (ELEV .GT. S) THEN 

*** EQ. (81 ) *** 

PINT=DSQRT ( 2 . *PI )*K1 *PRT6*PRT4A*PRT5B/ 

1 DSQRT ( PRT7A-PRT7B ) 

ELSE 
C 

C *** EQ. (82) *** 

C 

PINT=DSQRT ( 2 . *PI )*K1 *PRT6*PRT4B*PRT5A/ 

1 DSQRT (PRT7B-PRT7A) 

ENDIP 

PRINT 15, ' INTEGRAL FOR THE DIRECT WAVE IS', PINT 
15 FORMAT (A, 2P1 2. 6) 

GO TO 200 
C 

100 ARG3=(3.*LANDA/2.*G32C(ZER0,CMPSAD))**(2./3.) 

CALL H(ARG3,H1 ,H2,H1P,H2P) 

T1 =1 . -DCMPLX ( 0. 0, 1 . 0) /2 . * ( SPEED/OMBGA ) *ZI* 

1 CGPIZ2 ( ZERO , CMPSAD ) /CGPIZ ( ZERO , CMPSAD ) 
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T2=DCMPLX(0.0,1 .0)*(SPEED/0MEXJA)*ZI*(3.*LANDA/2. )** 

1 (2./3.)*CGPIZ(ZER0,CMPSAD) 

R=- ( T1 *H1 +T2*H1 P ) / ( T1 *H2+T2*H2P ) 

PRT8=€DEXP( PI/6. *DCMPLX(0. 0,-1 . ))/(CDEXP(DCMPLX(0.0, 1 . )*PI/6. )- 
1 DCMPLX(0.0.1 .0)*R) 

IP (IDEC .EQ. 3)G0 TO 60 
IP (IDEC .EQ. 2) THEN 

GO TO (30,40,30,40,50,50) REGION 
ENDIF 

IP (REGION .EQ. 5) GO TO 30 
IP (REGION .EQ. 6) GO TO 40 
C 

30 PRT7C=2.*G32PI2( ZERO, SADDLE) 

*** EQ. (83) *** 

PINT=DSQRT ( 2 . *PI ) *K1 *R*PRT6*PRT4B*PRT5B/ 

1 DSQRT ( PRT7C-PRT7B-PRT7A ) 

PRINT 15, ’ INTEGRAL FOR THE REFLECTED WAVE IS', PINT 
GO TO 200 

*** EQ. (85) *** 

0 PINT=DSQRT(2.*PI)*K1 *PRT6*PRT4B*PRT5B/DSQRT (-PRT7A- 

1 PRT7B)*PRT8 

PRINT 15, ' INTEGRAL FOR THE REFRACTED WAVE IS' .PINT 
GO TO 200 

50 PRT 6A=~CDEXP ( DCMPLX ( 0 . 0 , -1 . 0 ) *SADDLE*HDIST ) 

*** EQ. (84) *** 

PINT=DSQRT (2 . *PI )*K1 *PRT6A*PRT4B*PRT5B/DSQRT (PRT7A+ 

1 PRT7B)*PRT8 

PRINT 15, ' INTEGRAL FOR THE REFRACTED WAVE IS '.PINT 
GO TO 200 

60 PRT6B=CDEXP (DCMPLX (0.0, -1 . 0 ) * ( SADDLE*HDIST-PI /4 . ) ) 

PRT7C=G32PI3 ( S , SADDLE ) 

PRT7D=G32PI 3 ( ELEV , SADDLE ) 

FPI3=(2./(PRT7C+PRT7D) )**(1 ./3* ) 

CALL EVALB4 (S, ELEV, ANGLE, B0UND4) 

DELR=HDIST-B0UND4 

IF ((REGION .EQ. 6) .OR. (REGION .EQ. 4)) THEN 
ARG=-ABS(FPI3 # DELR) 

ELSE 

ARG=ABS ( FPI 3*DELR ) 

ENDIF 

CALL CGBAIR(ARG,AI,BI, AIP.BIP) 

*** EQ. (89) *** 

PINT =2 . *PI *K1 *PRT6B*PRT4B*PRT5B«FPI 3*AI *PRT8 



PRINT 15, ' INTEGRAL FOR THE WAVE BY CAUSTICS IS', PINT 


RETURN 

END 


SUBROUTINE EVALB4(S,ELEV,ANGLE4,BOUND4) 

*** THIS SUBROUTINE EVALUATES THE LOCATION 0^ *** 
***■ BOUNDARY IV *** 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON / ONE/ SPEED , OMEGA , ALPHA , DTOT 

GUESS2= ( OMEGA/ SPEED ) *DSQRT ( 1 ./(I .+DTOT)) 

IF (ELEV .LT. 1.E-5) THEN 

ANGLE4=DAC0S ( DSQRT ( ( 1 . +ALPHA*S+DTOT ) / ( ( 1 . +DTOT ) 
*(1 .+ALPHA*S)))) 

B0UND4=FNRAY ( ELEV , ANGLE4 )+FNRAY ( S , ANGLE4 ) 

GO TO 50 
END IF 

IF (ELEV .GT. S) THEN 

GUESS 1 =( OMEGA/SPEED )*DSQRT ( ( 1 .+ALPHA*S)/ 

( 1 . +ALPHA*S+DTOT ) )-0 . 0000001 
ELSE 

GUESS 1 = ( OMEGA/SPEED )*DSQRT ( ( 1 . +ALPHA*ELEV ) / 

( 1 . +ALPHA*ELEV +DTQT ) )- 0 . 0000001 
MT IF 
IC0UNT=0 
IC0U1JT=IC0UNT+1 
IF (ICOUNT .GT. 1 00)THEN 

PRINT*, ’SOMETHING WRONG WITH EVALB4! ! ! ' 

GO TO 50 
ENDIF 

Y1 =G32PI2 (ELEV , GUESS 1 )-fG32PI2(3,GUESS1 ) 

Y2=G32PI2 (ELEV , GUESS2 )-KJ32PI2 ( S , GUESS2 ) 

CALL ROOT ( GUESS 1 ,GUE3S2,Y1 ,Y2,BETA) 

IF (BETA .EQ. 0.0)G0 TO 10 

ANGLE4=DAC0S ( ( BETA*3PEED ) / ( OMEGA*DSQRT ( ( 1 . + 

ALPHA*S)/(1 . +ALPHA*S+DTOT ) ) ) ) 

B0UND4=FNRAY (ELF/ , ANGLE4 )+FNRAY(S, ANGLE4 ) 

RETURN 

END 




SUBROUTINE ROOT(XL,XR,FXL,FXR,SOLN) 
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*** THIS SUBROUTINE UTILIZES THE BISECTION *** 
*** METHOD FOR ROOT FINDING *** 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

S01fl=0.0 

IF ((FXL*FXR) .GT. 0.0)THEN 
XR=XL 
XL=TEMP 
XL' (XL+XR)/2.0 
GO TO 10 
ENDIF 

IF (ABS(XR-XL) .GT. 1.E-10)THEW 
TMP=XL 

XL=(XL+XR)/2.0 
GO TO 10 
ENDIF 

SOLN= (XLfXR ) /2 . 0 

10 RETURN 
END 


SUBROUTINE H(Z,H1 ,H2,H1 P,H2P) 

*** H USES SUBROUTINE CGBAIR TO CALCULATE 1/3 ORDER HANKLE *** 
*** FUNCTIONS FROM AIRY FUNCTIONS. *** 

IMPLICIT REAL*8 (A-H.O-Z) 

C0MPLEX*1 6 Z,AI,BI,AIP,BIP,K,KS,H1 ,H2,H1 P,H2P,ARG,CI 
CI= DCMPLX(O.DO, 1 .DO) 

PI= 3-H1592654DO 

ARG= DCMPLX(O.DO,-PI/6.DO) 

K= (12.D0)**(1 .D0/6.D0)*CDEXP(ARG) 

KS= DCONJG(K) 

CALL CGBAIR(-Z, AI,BI f AIP,BIP) 

HI = K*(AI-CI*BI) 

H2= KS*(AI+CI*BI) 

HI P= -K*(AIP-CI*BIP) 

H2P= -KS*(AIP4CI*BIP) 

RETURN 

EJD 

SUBROUTINE CGBAIR(Z, AI,BI, AIP.BIP) 

IMPLICIT REALMS (A-H,0-Z) 

CALCULATE AIRY FUNCTIONS FOR COMPLEX* 1 6 ARGUMENT 

REP. HANDBOOK OF MATHEMATICAL FUNCTIONS, ABRAMCWITZ AND STEGUN. 
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ENTRY 

CALCULATE ARGUMENT(Z) AND ABSOLUTE VALUE(Z) 

4 IP /Z/ UT 5 

THEN USE EQS. 10.4.2 THRU 10.4*5 FOR AI,BI,AIP,BIP 
10 EISE IP ARG(Z) LT Fl/3 

THIN CALCULATE ZETA(Z) 

USE EQS. 10.4.59, 10.4.61, 10.4.63, 10.4.66 FOR AI,BI,AIP,BIP 
20 ELSE CALCULATE ZETA(-Z) 

USE EQS. 10.4.60, 10.4.62, 10.4.64, 10.4.67 FOR AI,BI,AIP,BIP 
ENDIF 
MDIF 
EXIT 

end 


C0MPLEX*1 6 Z,AI,BI,AIP,BIP,ZETA,CZETA,Z14,SUM1 ,SUM2,SUM3,SUM4, 
1 ZETAP,FACT1 ,FACT2,SN,CS,FTERM,FPTERM,GTERM,GPTERM,F,FP,G,GP,Z3 
DIMENSION C (21 ),D(21 ) 

DATA Cl ,C2,PIRT,PI4/.355O23O539D0, .25881 94038D0.1 .772453851DO, 
+ .7855981 635DO/ 

+ .0371 35487654321 DO , . 0379950591 27800D0 , 

1 .0576491 9041 2669DO,. 11 609906402551 DO, 

+ . 291 591 59925074D0, . 87766696950998DO, 

2 5* 0794550501 751 IX), 12. 541 575552545DO, 

+ 55- 62273556591 4D0, 278. 46508077759D0, 

5 1555.1 6945201 27D0, 9207. 2065997258D0, 

+ 59892. 51 55658751)0,41 9524-8751 1655DO, 

4 51 48257. 41 78666DO, 251 9891 9-871 601 DO, 

+ 21 4288056. 96566DO,1 929575549. 1825D0, 

5 1 8355766957. 889D0/ 

DATA D/1 .DO, 

+ -.097222222222221 DO, -.0438850308641 97D0,-.042462830789894D0, 

1 -.0626621 63492031 DO,-. 1 241 0589602727D0,-. 508253764901 07D0, 

2 -. 92047999241 291 DO, -3. 21 04935846485D0,-1 2.8072930807551)0, 

3 -57. 50830351391 IDO, -287. 0332371 0920D0,-1 576.3573033370D0, 

4 -9446 . 3548230953DC , -61 335 . 706663847DO, -428952 . 40040004D0, 

5 -321 4536. 521 4006D0, -25697908. 383909D0,-21 8293420.8321 4D0, 

6 -1 963523788. 9909D0, -18643931 088. 105D0/ 

ABSZ=ABS(Z) 

IF(ABSZ.EQ.O) GO TO 3 

IF(ABS(DIMAG(Z))-LE.1 .D-12.AND.DREAL(Z).LT.0.D0) GO TO 5 
ARGZ=AT AN2 ( DIMAG ( Z ) , DREAL ( Z ) ) 

GO TO 4 

3 ARGZ=O.DO 
GO TO 4 

5 ARGZ=3 • 1 41 5926535898DO 

4 CONTINUE 

IF(ABSZ.GT.4.5DO) GO TO 10 

ASCENDING SERIES 
EQS. 10.4.2,10.4.3 


Z3=Z**3 
FTERM=1 .DO 
FPTERM=Z*Z/ 2 . DO 
GTERM=Z ' 

G?TERM=1 .DO 
GLIM=1 .D-13*ABSZ 
F=FTERM 
fp=fpt2w 

G=G7ERM 
GP=GPTERM 
DO 1 1=1,100 
13=3*1 

FTERM=FTERM*Z3/ ( ( 13-1 .D0)*I3) 

FPTIHM=FPTIRM*Z3 / ( 13* ( 13+2 . DO ) ) 

GTERM=GTIEM*Z3/(I3*(I3+1 .DO)) 

GPTE?M=GPTERM*Z3/ ( (13-2 . DO ) *13 ) 

F=F+FTERM 

FP-FP+FPTERM 

G=G+GTERM 

GP=G?4GPTEW 

IF(ABS(GTERM).LE.GLIM) GO TO 2 

1 CONTINUE 
PRINT 6000, Z 

6000 FORMAT ( / ' Z='2E1 4-5, ' ERROR IN CGBAIR, NONCONVIRGENCE ' ) 

2 AI=C1*F-C2*G 
AIP=C1 *FP-C2*GP 

BI=1 .732050808D0*(C1 *F+C2*G) 

BIP=1 . 732050808DO* (Cl *FP+C2*GP ) 

GO TO 9999 

C ASYMPTOTIC EXPANSIONS FOR /Z/ LARGE 

10 SIGN=1 .DO 
SUM1 =0.D0 
SUM2=0.D0 
SUM3=0.D0 
SUM4=0.D0 

C PI/3 = 1.047197551 

IF(ABS(ARGZ).GE.1 .3D0) GO TO 20 
C 

C /ARG(Z)/ LE PI/3 

C EOS. 10.4.59, 10.4.61, 10.4.63, 10.4.66 

C 

ZETA=CZETA ( ABSZ , ARGZ ) 

DO 11 1=1,12 
K=I-1 

ZETAP=ZETA**K 
SUM1 =SUM1 +SIGN*C ( I ) /ZETAP 
SUM2=SUM2+SIGN*D ( I ) /ZETAP 
SUM3=SUM3+C ( I ) /ZETAP 
SUM4=SUM4+D ( I ) / ZETAP 

1 1 SIGN=-SIGN 

Z 1 4=ABSZ** . 25DO*DCMPLX ( COS ( ARGZ /4 . DO ) , SIN ( ARGZ /4 . DO ) ) 
FACT1 = . 5D0*EXP (-ZFTA ) / ( PIRT*Z 1 4 ) 

FACT2= . 5D0*EXP (-ZETA )*Z1 4/PIRT 
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AI=FACT1 *SUM1 
A I P=-FACT2*SUM2 

* FACT 1 =EXP ( ZET A ) / ( PIRT*Z 1 4 ) 

FACT2-EXP ( ZETA ) *Z 1 4/PLRT 
BI=FACT1 *SUM3 
BIP=FACT2*SUM4 
GO TO 9999 

/ARG(Z)/ GT PI/3 

EQS. 10.4.60, 10.4.62, 10.4.64, 10.4.67 

20 ARGZ=ATAN2 ( -DIMAG ( Z ) , -DREAL ( Z ) ) 

ZETA=CZETA ( ABSZ , ARGZ ) 

DO 21 1=1 ,10 
K2=(I-1 )*2 
J=K2+1 

ZETAP=ZEIA**K2 
SUM1 =SUM1 +SIGN*C ( J ) /ZETAP 
SUM2=SUM2+SIGN*C ( J+1 )/(ZETAP*ZETA) 
SUM3=SUM3+S T ®*D(J)/ZETAP 
SUM4=SUM4+SIGN*D (J+1 )/(ZETAP*ZETA) 

21 SIGN=-SI(21 

Z 1 4= ABSZ** . 25DO*DCMPLX ( COS ( ARGZ/4 . DO ) , SIN ( ARGZ /4 . DO ) ) 
FACT1=1.D0/(PIRT*Z14) 

FACT2=Z1 4/PIRT 
SN=SIN(ZETA+PI4 ) 

CS=C0S(ZETA+PI4) 

AI=FACT1 *(SN*SUM1 -CS*SUM2 ) 

AIP=-FACT2* ( CS*SUM3+SN*SUM4 ) 

BI=FACT1 *(CS*3UM1 +SN*SUM2) 

BIP=FACT2* ( SN*SUM3-CS*SUM4 ) 

9999 RETURN 
END 


FUNCTION CZETA( ABSZ, ARGZ) 

IMPLICIT REAL*8 (A-H,0-Z) 

C0MPLEX*16 CZETA 
ARG=ARGZ*1 . 5D0 

CZETA= ( ABSZ** 1 . 5D0 ) *DCMPLX ( COS ( ARG ) , SIN ( ARG ) ) * . 66666666666667DO 

RETURN 

MD 


» 


FUNCTION FNRAY (Z, THETA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
COMMON /CNE/SPEED , OMEGA , ALPHA , DTOT/TWO/S 
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A=ALPHA*(1 . +ALPHA*S+DTOT-DCOS ( THETA ) **2* ( 1 .+ALPHA*S)) 
B=1 .+ALPHA*S+DT0T-DC0S( THETA )**2*(1 .+DT0T)*(1 .+ALPHA*S) 
IP (ABS(B) .LI. 1E-15)B=0.0 
C=DC0S ( THETA ) **2* ( 1 .+ALPHA*S)* ALPHA 
d=(i .+dtotWi .+ALPHA*3)*DC0S( THETA )**2 
E=C*(A*Z+B)/(A*(C*Z+D) ) 

P=DSQRT(C*(A*Z+B)/(A*(C*Z+D) ) ) 

PNRAY=DSQRT ( (A*Z+B)*(C*Z+D ) ) /A+ ( A*D-B*C ) *DL0G ( ( 1 +P ) / 
(1-P))/(A*DSQRT(A*C)*2.) 

RETURN 

END 


FUNCTION G32(Z,BETA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT 
AA=1 . +ALPHA*Z+DTOT 
BB=1 . - ( SFEED*EETA/ OMEGA ) **2 
ROOTBB=DSQRT ( BB ) 

PHI =DSQRT ( ( BB*AA-DTOT ) / ( BB*AA ) ) 

FOO=DLOG( (1 .+PHI)/(1 .-PHI)) 

G32=DSQRT (BB*AA**2-DT0T*AA )-0 . 5*DT0T/R00TBB*F00 

RETURN 

MD 


FUNCTION G32PI(Z,BETA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z^ 

COMMON /ONE/SPEED , OMEGA , ALPHA , DTOT 
AA=1 . +ALPHA*Z+DTOT 
EB=1 .-(SPEED*BETA/0MEGA)**2 
ROOTBB=DSQRT ( BB ) 

QUAN=EB*AA-DTOT 
IF (QUAN .LT. 1.E-15) THEN 
PHI =0.0 
CC=0.0 
ELSE 

PHI =DSQRT ( QUAN/ ( BB* AA ) ) 

CC=DSQRT (EB*AA**2-DT0T*AA ) 

ENDIF 

F00=DL0G((1 ,+PHI)/(1 .-PHI)) 

DD=0. 5*DT0T/R00TBB*F00 

G32PI=( ( -SPEED*BETA ) / ( ALPHA *OMEGA*BB ) )*(CC+DD) 

RETURx 

END 
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FUNCTION G32PI2(Z,BETA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /ONE/SPEED , OMEGA , ALPHA , DTOT 
AA=1 . +ALPHA*Z+DTOT 
BB=1 .-(SPEED*BETA/0MEGA)**2 
CC=1 . +2 . * ( SPEED*BETA / OMEGA ) **2 
ROOTBB=DSQRT ( BB ) 

PH1=DSQRT ( (BB*AA-DTOT )/(BB*AA ) ) 

F00=DL0G( (1 .+PHI)/(1 .-PHI)) 

G32PI2= (-SPEED/ ( ALPHA*0MEGA*BB**2 ) )* 

1 ((( BB* AA ** 2-CC *DTOT *AA ) /DSQRT ( BB*AA**2-DT0T*AA ) ) + 
1 ((CC*DT0T)/(2.*R00TEB))*F00) 

RETURN 

END 


FUNCTION G32FI3(Z,BETA) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT 

A= (SPEED*BETA / OMEGA ) **2 

AO=A/EETA 

A1 =1 .-A 

A2=1.+2.*A 

A3=5.-2.*A 

B=1 .+ALPHA*Z+DTOT 

P=DSQRT( (A1 *B-DTOT)/(A1 *B) ) 

F=DLOG((1 .+P)/(1 —P)) 

PRTO=- ( SPEED/ ( ALPHA*OMEGA ) )/A1 **2 

PRT1 =AO*B**2* ( 2 . *A*DTOT*B+4 . *DT0T**2-A 1 *B**2-3-*DTOT*B) 
PPT2=(A1 *B**2-DT0T*B)*"(3./2. ) 

PRT3=F*AO*A3 / A1 ** ( 3 • /2 . ) -2 . *A0*A2*DSQRT ( B ) / 

1 ( A1 *DSQRT ( A1 *B-DT0T ) ) 

PRT4=4.*A0/A1 *G32PI2(Z,BETA) 

G32PI3=PRTO*(PRT1 /FRT2-K) . 5*DT0T*PRT3 ) +PRT4 

RETURN 

END 


FUNCTION G32C(Z,BETA) 

IMPLICIT REAL*8 (A-H.O-Z) 

COMPLEX*! 6 BETA, PHJ, SORT 1 ,SQRT2,FOO,G32C,MODLOG 
COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT 
AA=1 .+ALPHA*Z+DTOT 

PHI=SQRT1 ( BETA, Z)/(SQRT2( BETA )*DSQRT(AA)) 
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F00=M0DL0G((1 . +PHI ) / ( 1 . -PHI ) ) 

G32C=SQRT(AA)*SQRT1 (BETA f Z)-.5*DTOT*POO/SQRT2(BETA) 

RETURN 

END 




FUNCTION SQRT1 (3ETA,Z) 

IMPLICIT REAL*8 (A-H.O-Z) 

C0MPLEX*1 6 BETA, AA, SORT 1 

COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT 

AA=1 — ( ( SPEED/OMEGA ) *BETA ) ** 2 . 

EB=1 .+ALPHA*Z+DTOT 

BRAN 1 = ( OMEGA/SPEED )*DSQRT ( ( 1 .+ALPHA*Z)/(1 . +ALPHA*Z+DTOT ) ) 
IF ((DREAL(EETA) .GT. BRAN1 ) .AND. 

1 (DIMAG(BETA) .LE. 0.0)) THEN 
SQRT1 =-CDSQRT ( AA*BB-DT0T ) 

ELSE 

SQRT1 =CDSQRT (AA*EB-DT0T ) 

END IF 
RETURI'I 
END 


FUNCTION S0RT2(BETA) 

IMPLICIT REALMS (A-H.O-Z) 

C0MPLEX*1 6 BETA.AA, SQRT2 

COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT 

A=ALPHA+DTOT 

AA=1 .-((SPEED/uMEGA)*BETA)**2. 
BRAN2=0MBGA/SPEED 
IF ((DREAL(BETA) .GT. BRAN 2) .AND. 
1 (DIMAG (BETA) .LE. 0.0)) THEN 
SQRT2=-CDSQRT ( AA ) 

ELSE 

SQRT2=CDSQRT ( AA ) 

END IF 
RETURN 
END 


FUNCTION MODLOG(QUAN) 

IMPLICIT REAL*8 (A-H,0-Z) 

C0MPLEX*16 QUAN.MODLOG 

IF ((DREAL(QUAN) .LT. 0.0) .AND. (DIMAG(QUAN) 


» 
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1 .GE. 0.0)) THEN 

MODLOG=LOG ( QUAN ) +DCMPLX ( 0 . 0 , -2 . «3. HI 5927 ) 
ETfTR 

MODLOG=LOG ( QUAN ) 

END IP 
RETURN 
END 
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FUNCTION GPIZ(Z,BETA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /ONE/SPEED , OMEGA , ALPHA , DTOT 
A=1 . - ( SPEED*EETA/ OMEGA ) **2 
B=1 . +ALPHA*Z+DTOT 
C=G32(Z,BETA)**(2./3.) 

GPIZ= ( 2 . *ALPHA/3 . ) *DSQRT ( ( A*B-DT0T ) / ( C*B ) ) 

RETURN 

END 


FUNCTION GPIZ2(Z,BETA) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /ONE/SFEED, OMEGA, ALPHA, DTOT 
A=1 . - ( SPEED*BETA/ OMEGA ) **2 
B=1 .+ALPHA*Z+DTOT 
C=G32(Z,BETA)**(2./3.) 

GPIZ2=(ALFHA**2./3. ) /DSQRT ( C* ( B**2* A-DT0T*3 ) )* 
1 ( DTOT/B- ( A*B-DTOT ) * (GPIZ ( Z , BETA ) / (ALPHA*C ) ) ) 
RETURN 
END 


ETJNCTION CGPIZ(Z,CBETA) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
COMMON /ONE/SPEED, OMEGA, ALPHA, DTOT 
COMPLEX* 1 6 CBETA , A , C , G32C , CGPIZ 
A=1 .-(SPEED*CBETA/0MEGA)**2 
B=1 . +ALPHA*Z+DT0T 
C=G32C ( Z , CBETA ) ** ( 2 . / 3 . ) 

IF (DIMAG(C) .EQ. 0.0) THEN 
SI=1 .0 
EI£E 
SI =-1.0 
ENDIF 


o o o o o 


CGPIZ=SI* ( 2 . * ALPHA/3 • ) *SQRT ( ( A*B-DTOT ) / ( C*B ) ) 

RETURN 

MD 
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FUNCTION CGPIZ2 ( Z , CBETA ) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /ONE/SPEED , OMEGA , ALPHA , DTOT 
COMPLEX*1 6 CBETA ,A f C, G32C , CGPIZ , CGPIZ2 
A=1 . - ( SPEED*CBETA/ OMEGA ) **2 
B=1 . +ALPHA*Z+DTOT 
C=G32C(Z,CBETA)**(2./3. ) 

CGPIZ 2= ( ALPHA**2 . /3 . ) /CDSQRT ( C* ( B**2*A-DT0T*B ) ) * 
1 (DTOT/B-(A*B-DTOT)*(CGPIZ(Z,CEETA)/(ALPHA*C) )) 
RETURN 
END 
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ABSTRACT 


In this study, the wave propagation over a finite impedance ground 
under a temperature lapse condition is analyzed. The solution is ex- 
pressed in terms of an integral. A surface wave term is found as a 
result of a pole occurring in the integration. The behaviors of the 
surface wave for different temperature gradients and ground impedance 
values are studied. This surface wave is seen to be important at low 
source frequencies. A grouna wave-like term is also seen in the total 
pressure field, which is the dominant term well inside the shadow re- 
gion. The calculated results of the pressure field show a very good 
agreement with the empirical model developed. 
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NOMENCLATURE 


a Sound speed 

a* Sound speed at z = ® 

T Temperature 

Tod Temperature at z = ® 

AT Temperature change between z = 0 and z = » 

G Function defined by (5) 

G Hankel transform of 3 

g Function define'' b*' (23) 

hj , h£ Modified Harkel functions 
i /-I 

P Pressure 

q A constant giving the source strength 

s Source height 

t Time 

Rg, Rj Reflection and/or refraction coefficient defined by (18) and 
(19) 

Tj , 1'2 Constants defined by (20) and (21) 

u) Circular frequency 

K, K 3 Constants defined (17) and (45) 

6 Hankel transform variable 

Z Ground impedance 

n Function defined by (22) 

z Receiver height 


Receiver radius 


Parameter defined 
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1. INTRODUCTION 


The prediction of acoustic noise is of increasing interest, since 
the control and reduction of noise in the community has become a major 
social and technical problem. It is particularly true for the area 
around an airport. An airplane which is taking off or landing is 
one of the most significant noise sources in the community. For over 
three decades, people have been trying to understand the physics of 
the generation, propagation, and pe.ception of aircraft noise. This 
thesis is concerned with only the propagation. 

The simplest physical model of outdoor sound propagation is a point 
source in free space which ignores the effect of ground, meteorolo- 
gical effects, and atmospheric absorption of sound. In the early 
1950s Ingard [1] and Lawhead and Rundick [2, 3] each obtained an analy- 
tical solution of the sound propagation from a point source in a homo- 
geneous atmosphere over a finite impedance ground. A finite impedance 
ground surface changes both the amplitude and phase of a sound wave 
being reflected from it. These solutions include the contributions 
of the direct and reflected waves plus a term now called the ground 
wave, and they describe an interference pattern above the ground sur- 
face. Although these solutions were accepted in this field little 
use was made of them for over twenty years. Because these models were 
really too simple, inconsistencies exist between them and the practical 
measurements made in laboratories and outdoors. 

In 1972 Wenzel [A] showed that a term, the so-called surface wave, 
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was missing in the previous solutions. After solving the same problem 
as Ingard, and Lawhead and Rundick by a different approach, he suggested 
that some of these inconsistencies could be the result of a surface 
wave which appeared to be an essential component of the total sound pres- 
sure field generated by a point source above a finite impedance ground. 
He further pointed out that under certain limiting conditions his 
results differed from those of Ingard by just this surface wave term. 
Chien and Soroka [5] also obtained a similar kind of asymptotic solu- 
tion for various limiting cases. However, their results applied to a 
more general case than that of Wenzel. They also found that in certain 
situations their solutions differed from Ingard by a surface wave 
term. This was also shown by Thomasson [6] and others. Tracing through 
the procedure used by Ingard, Thomasson found that Ingard, in the pro- 
cess of deforming the integration path while using the saddle-point 
method, failed to take proper account of a pole in the integrand. It 
is this pole which contributes the surface wave. 

In 1973, Embleton, Piercy and Olson [7, 8] introduced the term 
ground wave to acoustic propagation. They showed that the ground 
wave term was the significant term in the sound field near the ground 
and far from the source. Since the ground wave decays as the square 
of the horizontal distance from the source, it decays more rapidly 
than the direct and reflected waves but dominates at large distance 
because of the cancellation of the direct and reflected waves. But 
the physical interpretations of the surface wave and ground wave are 
still not clear. 

It should be noticed that all the solution models mentioned above 
are based on the assumption of homogeneous atmosphere. The refrac- 
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t i on effects due to the inhomogeneity of the atmosphere had not been 
considered. On the other hand, most of the available experimental 
data showed a difference in the sound pressure level between the homo- 
geneous and inhomogeneous atmosphere. For example, in 1959, Wiener and 
Keast [9] conducted a number of measurements of the difference between 
the simple source model and actual sound level, the so-called excess 
attenuation. They found that under the conditions producing an acoustic 
shadow zone the excess attenuation increased rapidly away from the 
shadow boundary and then leveled off to an approximately constant 
value well inside the shadow zone. This phenomenon could not be ex- 
plained by the models used then. Parkin and Scholes [10, 11], after 
a series of experiments in two different locations, reported that the 
ground excess attenuation could sometimes be negative, which differs 
from the theory even more. Pao and Evans [12] showed similar evidence 
of the disagreement between theory and experiment. More recently, 
Joriasson [13], on the basis of a similar study, concluded that the 
agreement L-tween theory used then and experiment was quite good for 
the case of soft boundary, but for the case of a harder boundary there 
was less agreement. 

Therefore, it is necessary to study the wave propagation behaviors 
in an inhomogeneous atmosphere, especially the behaviors of surface 
wave and ground wave under the conditions creating acoustic shadow 
zone. This is the purpose of the present project. 

The two major factors which may produce the shadow zone are the 
wind and temperature gradients. It is more difficult to find the solu- 
tion in the presence of the wind since the problem is not symmetric. 
But the effect of wind vector on the acoustic sound field is similar 


to that of temperature gradient (temperature lapse effect is similar 
to upwind propagation and temperature inversion effect similar to 
downwind propagation). Therefore, the only focus of this study is on 
the temperature gradient effect. 


2. FORMATION OF THE SHADOW ZONE 


A ray is defined as the line perpendicular to the wave front in 
the absence of wind. The idea of the ray is used here to examine 
the mechanism of shadow zone formation in the presence of a tempera- 
ture gradient. 

There are two kinds of temperature gradient conditions: inversion 
and lapse. For the case of a temperature inversion where the tempera- 
ture and consequently, the acoustic speed increases with height, the 
sound rays from a source above a finite impedance ground are refracted 
downwards (Figure la). This downwards refraction results in an in- 
crease of acoustic energy near the ground, and no shadow zone is created 
near the ground surface. Temperature inversion usually occurs during 
the night, especially a clear summer night. This will not be the sub- 
ject of the present study. 

During the daytime when a temperature lapse commonly exists, the 
temperature and speed of sound decrease with an increase of height. 
Thus, the sound rays from a point source above the ground are refracted 
upwards (Figure lb). As the result of this upward refraction, there 
exists a ray which just grazes the ground at a distance Rg from the 
source. Beyond that distance, an acoustic shadow zone, into which no 
ray can penetrate, is observed. From Figure 1, the wave pattern in 
the presence of a temperature gradient, is seen to be symmetrical 
around the sound source. This symmetry is also the reason why the 
three-dimensional wave propagation problem can be simplified as a two- 


(a) i INVERSION 



Figure 1. Sound rays in the presence of temperature, 
(a) Inversion (b) Lapse 
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dimensional wave propagation problem. 

A mathematical model of the ray pattern in the presence of a tem- 
perature lapse has been established in the previous study [14] by 
applying Snell's law, which shows that the physical space can be divi- 
ded into seven regions (Figure 2). Each region is characterized by a 
different two-ray combination of direct, reflected and/or refracted 
rays, except for region seven which is the shadow region. It is also 
noticed that there are two kinds of refracted rays. To a receiver 

located at point A of Figure 3, refracted ray B has contacted caustic 
and ray C has not (Figure 3). 

In this model, a temperature profile: 

AT , % 

T = Too + ; (1) 

is assumed, where AT > 0 corresponds to lapse condition. This profile 
has been shown to fit the experimental data for a lapse condition very 
well [15]. 



Figure 2, Ray diagram under a temperature lapse condition. 



HORIZONTAL DISTANCE 




3. FORMAL SOLUTION 


The problem of a point source over a finite impedance ground in 
the presence of a temperature lapse can be stated as follows. We 
are looking for a function P(z,r,t) which satisfies the governing 
wave equation: 


1 _ 

a 2 (z) 



V 2 P = Q 


( 2 ) 


where Q represents the source function: 


Q = — 6(r) <5{z - s) e iljt 

ITT 


(3) 


and a(z) is sound speed profile: 


a(z) 



(4) 


which follows from the temperature profile (1). Ine solution must 
satisfy the boundary condition: 

P = - W Z (5) 

at z = 0, where Z is the ground impedance, W is the vertical component 
of the acoustic velocity. A "radiation condition," which requires the 
waves to outgoing as either r or z tends to infinity, is also used. 
Due to the symmetry described in the previous section, this will be a 
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two-dimensional wave propagation problem, z represents elevation of 
the receiver and r is the horizontal distance from the source. 
A solution developed by Van Moorhem [16] is used here. It gives: 

P = e iujt G(r,z) (6) 

where G is the function of r and z only and is defined as the inverse 

Hankel transformation of a function G, i . e . : 

00 

G = f G(r,z) r J 0 ( Br ) dr (7) 

C 

oo 

G = f G(B,z) 8 J 0 (Br) dB . (8) 

J 0 

G function in the complex plane is given as: 

G = K h 2 (n(B.z)) [h ( n( b,s ) ) + R 0 f >2 (n(B.s))] (9) 

in region A of the B plane (Figure 4a) 

G = K h 2 (n( 6 ,z ) ) Chj ( n( B.s ) ) + Ri h 2 ((n(B,s))] (10) 

in region B 

i 2 tt i 2 it 

G = K h 2 (n(B,z)) [hjfnU.sJe"^"" ) + h 2 ( n( S,s )e^~ )] (11) 

in region C 

i 2 TT i 2 TT i 2 TT 

3 = K h 2 (n(B,z)e 3 ) [hj(n(B,s)e 3 ) * R 1 h 2 ( n ( B.sje”'*”)] (12) 

in region D for z > s and 


G = K h 2 (n(B.s)) [hj(n(B,z)) + Ro h 2 (n(B,z))] 


(13) 


in region E (Figure 4b) 
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G = K h 2 (n(B'S)) [h!( n(s,z)) + Rj h 2 (n(B,z))] (14) 


in region F 


i 2 t i 2 tt 

G = K h 2 (n(B,s) ) [hj(n(8,z)e 3 ) + Rj h 2 (n(B,z)e 3 )] 


(15) 


in region G 

i 2* i 2 ti i 2 it 

G = K h 2 (n(B.s)e 3 ) [hj(n(B,z)e 3 ) + Rj h 2 (n(B,z)e 3 (16) 

in region H for z < s, where 


_Q 1 1_ 

!2i * 2/3 Vg 2 (S,s) Jg z [B,z) 

T 1 hi (n(B.O)) + i T 2 hj ( n( B.O) ) 
Tl h 2 (n(B.O)) + i T 2 h£ (n(B,0)) 

TT 

e ‘^6 

e^6 ~ TR 0 

1 Z 9 zz ( B,0 ) 

a X - - — 

2 pa. g z (B.O) 


(17) 


(18) 


(19) 


( 20 ) 


T 2 = — (2/3\) ?/3 g 2 (6,0) 

p 3 to 


( 21 ) 


n = ( 2/3 A ) 3/2 g(M) 


( 21 ) 
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The g function is defined by 


J 3/2 (B,z) = rj( 1 + az + ^) ((1 - % B 2 ) ( 

Too <i> 


at, at, 

1 + aZ + — ) - — ) 

T T 


1 £l 

2 To, 


'X/ 


8 

a) 


in ( 


1 ♦ ♦ . 
1 - t 


(23) 



The function g(z,s) has two branch points, one at 


U) 

B = B z = — 


'1 + aZ 


AT 

1 + aZ + — 

T„ 


(25) 


and the other at B = Bw = o»/a «• It is noticed that equations (11) to 
(18) contain the functions g(0,s), g ( s , B ) and g(z,B). Thus these solu- 
tions involve four branch points B 0 » 8$> 8z» and B^. These branch 
branch points form the boundaries of regions of different forms of the 
solution on the real B axis. The regions of validity of these solutions 
in the complex B plane are separated by the branch cuts of the different 
g functions. Each g function has two branch cuts. The first branch 
cut results from a square root in the g function and is chosen such 
that / -1 = -i, which is required to produce the exponentially decay 
behavior of the pressure in the far field. The second branch line is 
the result of taking the 2/3 power to calculate g from g 3/2 . 


It is also noticed that every form of the solutions consists of 
two terms. The first term involves no ground impedance and represents 
the direct wave term. The other term which includes the effect of the 
ground is the reflected or refracted wave term. Physically the four 
branch points 6 0 » 8 S , B z , and B^ represent the rays which have the turn- 
ing points at height of zero, s, z and infinity, correspondingly. Thus, 
the second term in the solutions represents a ray which is reflected 
by the ground for real B , where 0 < 8 < B 0 , a refracted ray which has the 
turning point above the ground and below the altitude of the receiver 
for real B and B 0 < 6 < 8 Z or a ray which has the turning point 
above the receiver and less than infinty for real B and B z < B < 8^. 
Those rays which have the turning points above the source, or real 8 
and Bs < 8 < By, are not seen in the ray diagram (Figure 3). For 
8 > 8^ or Im(8) * 0, the physical interpretation is not clear. 

These solutions are very complicated to evaluate since the inverse 
Hankel transformati on has to be carried out. Some approximating methods 
have to be adopted. The approximation used here is the saddle point 


method. 


4. SADDLE POINT METHOD 


The saddle point method has been frequently and successfully used 
in research on the wave propagation problems. It is an approximate 
integration method for integrals of the form [17, 18]: 

r vf(B) 

I(v) = J F(B)e d& (26) 
c 

where v is a large parameter and c is a contour in the complexe b plane. 
From complex variable theory, it is known that the integration path 
can be deformed while taking the proper account of the poles of the 
integrand. What then would be desirable is to deform the path so that 
there is very little contribution to the value of the integral over 
the whole path except in the immediate neighborhood of one or mere 
points on the path, where practically the entire contribution takes 
place. These points are found to be the saddle points, where 


9f 
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(27) 


Around these points the function surface looks like a saddle or a moun- 
tain pass. If the path is chosen such that the exponential term in (26) 
decreases rapidly as one moves away from the saddle points, the objective 
is achieved. This path is call a steepest descent integration path and 
is defined by requiring 
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Im ( f ( B ) ) = const. = Im (f(B sp )) (28) 

to eliminate oscillations in the integral and allow only decay. This 
results in the approximate integration. 


F(B 


sp>e Vf 'V 


+ ie 


(29) 


where 9 is the angle at which the path of the integration crosses the 
saddle point. If no saddle point has been found, the integral can 
be approximated as zero. 

Applying the saddle point approach to the present problem, the in- 
verse Hankel transform is first rewritten into the form 

oo 

6 ( r,z ) j G(B,z) B (fir) dB . (30) 

ao 

The exponential terms are obtained by assuming the argument of the Hankel 
functions are sufficiently large and the use of asympototic forms are 
valid. It then follows that for every point in the physical space, 
two saddle points occur, each corresponds to a specific ray which 
passes through that point, a direct, reflected or refracted ray. Extra 
caution has to be taken while tracing the integration path, since sev- 
eral branch cuts exist in the B plane. It will be possible for the 
integration path to run into a branch cut, which will require a term 
beyond the simple saddle point result. 

A detailed examination of the integration path has been made. For 
the direct wave (Region 1, 2 , 3, 4 of Figure 2), the integration path 
crosses the real s axis at the saddle point at an angle of tt/4 and the 
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whole path will not interfere with any branch cuts (Figure 5). For a 
reflected wave (Region 1, 3, 5), the path of integration also crosses 
the real B axis at the saddle point at an angle of tt/ 4 . The path is 
also simple for most cases, but under certain conditions the path 
may run into the branch cut of the g(0,B) function in the upper b plane, 
C] in Figure 6 , and moves along this branch line, then reappears on 
the same side of the branch line (Figure 6 ). Note that Van Moorhem [16] 
has shown that although ci is a branch line for g function, the inte- 
gral is continuous across it. For the refracted wave which has already 
contacted the caustic (Region 2, 4, 6 ), the path will run into both 
branch cuts of the g ( 0 , B ) function, ci and c? in Figure 7. The path 
starts at -®, and runs into branch line C 2 first and reappears right 
at the other side of this branch cut. It then passes through the saddle 
point at an angle of tt/4 and runs into the branch cut cj. Two possi- 
bilities exist after this: it might move along the cl line and reappear 
at the other side of cl (Figure 7): or it might move along the ci line 
and run into the branch point s^. In the latter case, a modified path 
has to be chosen such that the path will not run into the branch point 
B^ where the g function is infinite. This modified path is shown in 
Figure 8 and requires a constant phase of tt/4 as the path moves away 
from the point b u . 

For the refracted wave which has not contacted the caustic (Region 
5, 6 ), the path of integration is quite different. This path will cross 
the saddle point at an angle of 3 tt/ 4 rather than n/4 as in the other 
cases. The upper part of the path will turn right and merge into branch 
cut cl. The two situations discussed above then hold here. The lower 
part of the path will disappear into the branch line C 4 (Figure 9). 
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Figure 5. Integration path associated with direct and reflected 
waves. 









have contacted caustics. 
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► 


But. it is noticed that in Regions 5 and 6 , both saddle points result 
► from the second term of the solution. The integral for the first 

term, the direct wave term, will be approximated as zero since no 
saddle point has been found for this terni. Furthermore, since the 
two saddle points come from the same term, a single path has to be 
found such that coincides with most parts of the integration paths 
of the two saddle points and passes through these two points. The 
simplest path which has all these properties is shown in Figure 10. 
It consists of four parts: p^ coincides with part of the integration 

path of first saddle point: p£ is a straight line which connects both 
saddle points: P 3 coincides with part of the integration path of the 
second saddle point; and P 4 is the modified steepest descent path. The 
integration along path p£ can be done numerically. Physically, this 
integral represents the contribution of those rays radiated from an 
infinity height between the two rays represented by the two saddle 
pointr Since these rays do not pass through the source, these contri- 
butions are practically small and are ignored in this study. The saddle 
point approximation still holds for the integration along p^ and P 3 . 
The integration along the ^4 is done by expanding the integrand around 
point Bu and the result will be presented later in this section. 

For a point in the shadow region, we can still find two saddle 
points. These saddle points are complex and represent rays with complex 
angles. Their physical meaning is not clear, but the path of integra- 
tion for these saddle points have the same characteristics of the 
cases described above (Figure 11). 

Thus, we come to the final solutions. For direct waves with z > s. 


1 


Regions 1 and 2 of Figure 2, we get: 
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JTl T K3 

= / 2 
a 4> t 


J- 


36 * 


Tl 

e S P 2 h 2 h l (n(B S p.s)) (31) 


where 


♦ = _x g2/3 ( 3 ,s) + X g 2 / 3 (B.z) - - 


(32) 


and B sp is the root of 


3*, 

r +— 1 = 0 • 

38 


(33) 


For direct waves with z < s, Region 3 and 4 of Figure 2, we get: 


°D 


J2v K 3 -i(B r - -) 

L S P 2 h-. ( n( 8 .z)) h 2 ( n( 3 »s ) ) 

1 sp c sp 


2 

9 4>? 
3B 2 


(34) 


where 


4 » 2 = X g 3/2 (8 ,s) - X g 3/2 (6, z) -- ( 35 ) 

and B sp is the root of 

r + —2 = 0 . (36) 

38 

For reflected waves. Regions 1, 3 and 5 of Figure 2, we get: 
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J2v K 3 



■i ( Bsp r ~ 2^ 


(h 2 (n(8 sp ,z))h 2 (n(B sp ,s))) 


(37) 


where 

*3 = x g 3/2 ( 6 . s ) + X g 3/2 (B,z) -2 X g 3/2 (0,0) - y 
and B sp is the root of 


(38) 


r + 2*3. = 0 • ( 39 > 
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For refracted waves which have not contacted the caustic. Regions 4, 5 
and 6 of Figure 2, we obtain 


G 


JT* k 3 Rj 



V 


-i ( 6 r - it) 

e S P h 2 ( n( B s , z ) ) h ^ (n(B sp .s)) 

(40) 


where 


$ 4 = X g 3/2 (B,s) + X ( 0 » z ) + — tt 


3/; 


11 


(41) 


and 6 


sp 


r + 


is the root of 



90 


(42) 


For refracted waves which have contacted the caustic. Regions 2, 4 and 
6 of Figure 2, we get: 
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g R = 


VJwJCaJ Ll ^1% r - 2» h2 (n(B sp ,z)) h 2 (n(e sp .s)) 


r. 3^3 

36 2 


(43) 


where 


$5 = *4' 


(43) 


and B S p 
of (42). And 


is also a root (and the smaller of the two roots for Region 6) 


*3 s 


12 a 

— sp 


241 X 2 ' 3 J 9 Z (B sp .O J*Z 
If the modified path p4 is used, an extra term 
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(46) 
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has to be included. This term decays as -3/2 power of the horizontal 
distance but it will be the dominant term far into the shadow. Physi- 
cally, this term represents a g> ound wave-like term but it is not 
exactly a ground wave since the ground wave decays as -2 power of the 
horizontal distance. 

One problem with this approach is that at the shadow boundary we 

have 


ZB 98 


(47) 


and a singularity occurs in (40) and (43). To avoid this singularity, 
a method presented by Saches and Silbiger [19] is used, which gives a 
nonuniform asymptotic solution valid around the shadow boundary. This 
approach gives 


Wj 1/3 

G R = 2, X 3 (8 c ) R X (B c ) (^f 4 ) 


-i(B c r - 4 ) 


1/3 

Ai (( — j 4 ) Ar) h 2 (n(B c .z)) • h 2 (n(S c ,s)) 

9 8 


(48) 


where 8 c is the root of 


a 2 


98 


2 (+4 (B c ) ) = 0 


(49) 


and 


, 
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Ar = r - r c (50) 

where r c is the caustic radius 

a 

r c = - — (414 (Sq)) . (51) 

3B 

Reference 20 shows this solution is equivalent to the one-term series 
expansion of the uniform asymptotic formulation which can be approxi- 
mately used in the shadow region. 


5. SURFATE WAVE 


A surface wave is a wave whose amplitude exponentially decays with 
the height above the ground and decreases with the distance from the 
source as the product of a -1/2 power and an exponential. Thus, for 
a certain range of distance r, if the exponential term is weak, this 
i — 1/2 decrease of the amplitude will make the surface wave dominant 
in the pressure field near the qrcund. 

As suggested by most of the previous researchers, the surface 
wave is the result of the presence of the ground with complex impe- 
dance, and mathematically it results from the contribution of a pole 
which is located inside the path of integration. 

Examining the solution forms, Equations (9) through (16), it is 
found that only two terms involve the ground impedance. R 0 and Ri . 
An effort has been made to find the poles of these two terms. 

The method used here is called the modified Newton metnod which 
is valid to find the root of an equation of a complex variable 

F(0) = 0 (52) 

where B is a complex variable. Since the terms R 0 and R^ involve the 
ground impedance, either experimental data or a model must be used. 
The impedance model suggested by Chessell [21] is adapted here, which 
gi ves 


2 = R + i X 


(53) 


where 


R f -75 

* 1 + 9.08 (-) (54) 

P 3<jo 0 

1 f -.73 

__ = - 11.9 (-) (55) 

pa* o 

where a is the flow resistance and f is the frequency. This model has 
been shown to fit the experimental data very well. For the grass 
covered ground, such as one around an airport, the flow resistance 
a is suggested to be between 150-300 c.g.s. unit by Piercy and Embleton 
[ 22 ] . 

Numerical calculation shows that the poles of Rj are located inside 
the region A or E of the complex e plane in Figure 4 where the solu- 
tion forms involving Rj are not valid. Thus, these poles do not give 
any contribution to the total solution. 

On the other hand, the poles of R 0 are of great interest. It is 
initially noticed that when temperature gradient is zero, an exact 
solution of the wave Equation (2) is obtainable. A set of poles of 
this exact solution for different ratios of o/uj is then calculated, and 
so is a similar set of poles of R 0 for aT/T* near zero. A comparison 
shows that these two sets of poles are located very close to each 
other as would be expected. Thus, it appears that these are related 
and it is the surface wave producing poles. 

Pole of P 0 a } also found for larger values of AT/T*,. These poles 
are located at the right hand side of the branch cut C 3 (Figure 12). For 
certain combinations of AT/T,,,, a, a and u>, the poles are found to be 
located inside the integration path. Therefore the residues of the poles 
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have to be added into the total solution according to the Residue 
Theorem. The residue is given as 

G s = -U KB p H q ( 2 ) ( B p r ) h 2 (n(B p ,z)) h 2 ( n( B p ,s ) ) 

- Tj hi (n(fJ p .0) ) + i T 2 h i (n(B p .O)) 

. (56) 

— (Tj h 2 (n( B p ,0) ) + i T 2 h^ ( n( B p .O) ) 

or 1 using the asymptotic form 

TT 

f2Q~ -i Bn r + i — 

G s = - i v k/_P h 2 (n(B p ,z)) h 2 (n(B p ,s)) e 4 

- Tj hj (n( B p .0 ) ) + i T 2 hj (n(B p .O)) 

. (57) 

~ ( T 1 h 2 (n(B p ,0)) + i T 2 h 2 ( n( B p ,0))) . 

This residue is found to describe a surface wave which exponen- 
tially decays with the height from the surface (Figure 13) and is the 
product of an exponentially decreasing term, since the pole locates at 
the lower s plane, and the reciprocal of the square root of distance 
from the source. It is also noticed that these poles which might be lo- 
cated inside the integration paths are very close to the real 6 axis, 
i .e. , 


0 < -Im(B p ) << 1 . (58) 

Therefore, in the midrange of the horizontal distance, the exponential 
decay is quite weak which makes the surface wave a dominant component 
in the pressure field near the ground. But in the far field, the ex- 
ponential decay behavior will quench the surface wave (Figu'e 14). 
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The effects of the temperature profile and the ground impedance 
on the surface wave are also analyzed. For a fixed a, o and u, the lo- 
cation of the pole will move closer to the branch cut C 3 with an in- 
crease of AT/Too, thus providing less chance to produce the surface 
wave. For a fixed aT/T*, 0 and w, the chance to see the surface wave 
increases with the value of a. These two facts show that the surface 
wave is more likely to be seen in the homogeneous atmosphere than in 
the inhomogeneous atmosphere. On the other hand, for a given temperature 
profile, the chance to tee the surface wave and the magnitude of the 
surface wave increase with the ratio of o/m, which generally meets the 
criterion for the existence of a surface wave developed by Chien and 
Soraka in the homogeneous case, since the magnitude of the imaginary 
part of the ground impedance increases with the ratio of o/u> faster than 
the real part of the ground impedance. Therefore, for a given uj, the 
harder the surface, the stronger the surface wave and the more the 
chance to see the surface wave: or, for a given flow resistance, the 
higher the angular frequency, the weaker the surface wave and the less 
the chance to see the surface wave. Figure 15 is a plot of surface 
wave versus frequency at flow resistance of 300 c.g.s. unit, which 
shows the magnitude of the surface wave decreases with an increase of 
frequency very fast and it becomes negligible when frequency > 500Hz. 
This result is similar to that of Piercy, Emblet.on and Sutherland 
(Figure 11 of reference 23). 




6. RESULT OF CALCULATION 


In this section, some results of calculation are presented and 
the flow tesistance is assumed to be 300 c.g.s. units. First, it is 
noticed that the ray tracing model is based on the assumption of high 
frequency. But how high the frequency must be, so that this solution 
will be valid, is still a question. Figure 16 is a plot of sound 
pressure levc’ versus frequency at z = 4m, r = 20m so that the 
receiver is ’ocard in region 3 of Figure 2, which shows that the 
sound pressure level oscillates about a mean level of 46 dB. From 
this plot, it is noticed that no strange behavior occurs at very low 
frequency of about 50 Hz. 

Figure 17 is a similar plot but at r = 40m, it is noticed that 
a 6 db decay of the mean sound pressure level for doubling the distance 
is observed. The inteference pattern has a longer cycle than the pre- 
vious plot, and no strange behavior is obser.-ed at low frequencies, 
too. 

Figure 18 is another plot of sound pressure level versus frequency. 
The location of the receiver is chosen such that both the surface wave 
and the ground wave-like term are observed. This plot shows that the 
surface wave dominates the pressure field at low frequency. For fre- 
quency higher than about 5°0Hz, the surface wave has no effect on the 
sound pressure level at all. On the other hand, the ground wave-like 
term has less effect on the total sound pressure field at very low fre- 
quencies than the surface wave, but it shows a strong effect at 're- 



SOUND PRESSURE LEVEL VERSUS FREQUENCY 



Fiqure 16. Sound pressure level versus frequency with A T/T 
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Figure 17. Sound pressure level versus frequency with aT/T, 
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Fiqure 18. Sound pressure level versus frequency with AT/T 
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quencies higher than 500 Hz and less than 3000 Hz. 

Figure 19 is a plot of the ground wave-like term versus frequency. 
It shows that this term also decreases with the increase of the fre- 
quency, but a comparison between Figure 15 and Figure 19 indicates 
that the ground wave-like term decreases much slower with the increase 
of the frequency than ti.e surface wave term. 

Figure 20 is a plot of sound pressure level versus horizontal 
distance from the source. The regions in the physical space and the 
empirical model of Wiener and Keast [9] are also shown. The simple 
source model has been suggested by Wiener and Keast for the region 
between the source and the shadow boundary. Beyond the shadow bound- 
ary, a rapid increasing excess attentuation is applied to the simple 
source model and finally a large constant excess attentuation is 
applied to the simple source model when the receiver is well inside 
the shadow region. This empirical model is based on a number of experi- 
ments and summarizes the best experimental data collection available 
now. This plot shows a very good agreement between the present mathe- 
matical model and the empirical model. Since the simple source model 
only gives the mean sound pressure level, the oscillating behavior of 
the present model around the mean value is expected. 

Figure 21 is a more detailed plot of the previous one which shows 
how each component of the pressure field add together to yield the total 
sound level. It is noticed that inside the shadow region there are 
three terms contributing to the total sound pressure level: the surface 
wave term; the nonuniform asymptotic term; and the ground wave-like 
term. The surface wave term has almost no effect on the total sound 
pressure field since it has decayed in an exponential manner and the 
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frequency of the source is relatively high. The nonuniform asymptotic 
term which physically represents the energy diffracted from those rays 
which pass close to the shadow boundary is also small when the receiver 
is well inside the shadow region because the exponential decay of the 
Airy function. Equation (48), inside the shadow region. The dominant 
term which gives the saturated behavior of the excess attenuation 
is the ground wave-iike term. Again, this term is not exactly a ground 
wave term since it decays with the horizontal distance as a -3/2 power. 
The overshooting navior of the sound pressure level just inside the 
shadow region results from the interference of the three terms and is 
actually observed, sometimes, in the practical measurements (Figure 2 
of reference 9) . 

Figure 22 is also a detailed plot of the sound pressure level 
versus horizontal distance. It is observed that the surface wave has 
a strong effect on the pressure field for horizontal distance between 
30m to 600m. Beyond that distance the exponential behavior actually 
makes the surface wave negligible. 

Figure 23 is another plot of sound level versus horizontal dis- 
tance. The heights of the source and the receiver are the same as in 
the experiments done by Wiener and Keast. The saturated value of 30 
dB for excess attentuation far into the shadow suggested by them is 
observed. Figures 2<* through 26 are also the plots of sound pressure 
level versus horizontal distance at different heights and under dif- 
ferent temperature profiles, which show the similar kind of behavior 
to the previous plots. The saturated values of the e:.cess attenuation 
are found to vary with the height of the receiver, the temperature pro- 
file and the source frequency. The dependancy of the ground wave-like 
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Figure 24. Sound pressure level versus horizontal distance with 
AT/T. = .01, a = 6, io = 6,000, s = 3.67m, z = 1.52m. 
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Figure 26. Sound pressure It el versus horizontal distance with 
AT/Tc = 0.1, a = 6, u = 4,000, s = 3.67m, z = 4.67m. 
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term on the frequency has been shown in Figure 19. Figure 27 is a plot 
of the ground wave-like term versus height of the receiver. It is no- 
ticed that the pressure level of the ground wave-like term increases 
with the height of the receiver for height less than lm, it then de- 
creases with the continuous increase of the height. It also decays 
slower than the surface wave with the increase of the height (Figure 
14 and Figure 27). 

Since Wiener and Keast only conducted the experiments at a fixed 
receiver height and more than ±10 dB spread of the measured data exist 
in their experiments, this present mathematical model is believed to 
fit the experiments very well. 

Embleton, Piercy and Olson [7] suggested that the ground wave term 
will be more dominant in the pressure field if the source and the re- 
ceiver are both located very close to the ground. Figure 28, where the 
heights of the source and receiver are both .5 meters, shows a similar 
kind of behavior of the ground wave-like term. Over the whole hori- 
zontal distance range, this ground wave-like term dominate the sound 
pressure level, except for horizontal distance less than 100m, where the 
surface "he has a very strong effect. 







7. CONCLUSIONS 


In this project, the behavior of sound propagation under a tempera- 
ture lase condition has been studied, especially the behavior inside 
the shadow region. 

The surface wave has been found under certain combinations of 
temperature profile and ground impedance. Mathematically, it is the 
result of the pole which is located between the original integration 
path and the steepest descent integration path. This surface wave has 
been seen to be important in a range of the horizontal distance from 
the source when receiver is close to ground and the frequency is low. 
At large distance, no matter how low the frequency is the exponen- 
tial decay behavior of this surface wave will be strong enough to 
quench it. The surface wave might also be the reason why we sense 
the negative excess attentuation in a certain range of horizontal 
distance at low frequency, which is observed in the experiments. Also 
this study shows that the surface wave is more important in the homo- 
geneous atmosphere. 

A term which decays as -3/2 power of the horizontal distance is 
also found. This term is required mathematically to keep the steepest 
descent integration path from running into a branch point. It is 
this term that gives the saturated excess attenuation inside the 
shadow region which was observed by Wiener and Keast. This term is 
believed to be a ground wave-like term. 

The present mathematical model is seen not to give a strange 
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behavior at low frequency. More experiments are needed to prove the 
validity of this model at very low frequency. These can be achieved 
by carrying out the outdoor measurements over a wide range of fre- 
quencies, from several hundred to several thousand Hertz, under dif- 
ferent temperature profiles. 

A comparison between present model and the empirical model sug- 
gested by Wiener and Keast has been made, which shows a very good 
agreement between the theory and the experiments. 
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PROGRAM PRESSURE 




* 

THIS 

IS THE COMPUTER PROGRAM OF A MATHEMATICAL 

# 

* 

MODEL FOR EVALUATING THE SOUND PRESSURE LEVEL AT ANY 

* 

* 

POINT IN THE PHYSICAL SPACE UNDER A TEMPERATURE 

* 

* 

LAPSE CONDITION. THE PROGRAM CONSISTS OF 1 MAIN ROU- 

* 

* 

TINE, 18 SUBROUTINES AND 10 FUNCTIONS. 

* 

* 

FOLLOWING IS A LIST OF STANDARD SYMBLES USED IN 

* 

* 

* 

THIS PROGRAM: 

« 

# 

* 

S 

ELEVATION OF THE SOUND SOURCE 

* 

* 

Z 

ELEVATION OF THE RECEIVER 

* 

* 

HDIS 

HORIZONTAL DISTANCE OF THE RECEIVER FROM 

* 

* 


THE SOURCE 

* 

* 

Z1 

COMPLEX GROUND IMPEDANCE 

* 

* 

R1 

FLOW RESISTANCE 

* 

* 

OMEGA 

ANGULAR FREQUENCY OF THE SOURCE 

* 

* 

DTOT — 

AT/Too 

* 

* 

speed 

SPEED OF THE SOUND AT INFINITY 

* 

* 

ALPHA 

CONSTANT 

* 

* 

REGION 

THE REGION OF THE LOCATION OF THE RECEIVER 

* 

* 


IN THE PHYSICAL SPACE 

* 

* 

BRW 

BRANCH POINT Bw 

* 

* 

BRS 

BRANCH POINT Bs 

* 

* 

brz 

BRANCH POINT Bz 

* 

* 

BRO 

BRANCH POINT Bo 

* 

* 

SADPOT — 

SADDLE POINT ARRAY 

* 

* 

P 

SOUND PRESSURE STRENGTH 

* 

•* 

SPL 

SOUND PRESSURE LEVEL PREDICTED BY SIMPLE 

* 

* 


SOURCE MODEL 

* 

* 



# 




HHHHHf* 


* 

* 

* 

* 

# 

* 

* 

* 

* 

* 

* 

* 

* 

* 


LIST OF INPUT: 
DTOT: 
ALPHA: 
S: 
ZL: 
ZT: 
NZ: 


OMEGAL: 

OMEGAT: 

NO: 


AT /Too 

TEMPERATURE PROFILE CONSTANT 
HEIGHT OF SOURCE 
LOWER LIMIT OF RECEIVER'S HEIGHT 
HIGHER LIMIT OF RECEIVER'S HEIGHT 
NUMBER OF CALCULATING POINTS IN 
Z DIRECTION 

LOWER LIMIT OF SOURCE ANGULAR 
FREQUENCY 

HIGHER LIMIT OF SOURCE ANGULAR 
FREQUENCY 

NU MBER OF CALCULATING POINTS OF 
DIFFERENT FREQUENCIES 


* 


* 


ooooooooo 


61 


HDISL: CLOSEST DISTANCE OP RECEIVERS 
LOCATION 

HDIST: FARTHEST DISTANCE OP RECEIVERS 
LOCATION 

NH: NUMBER OP CALCULATING POINTS 
IN HORIZONTAL DIRECTION 


I I n mmxmmmmm m mmmmiimmmmmm 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

CGi PLEX*16 BETA,Z1 ,P1 ,P2,P3,P4,P5,P,G32PI,RC2 
COMMON /BRANCH /BRO , BRS , BRZ , BRW 
COMMON /CONSTANT/SPEED , OMEGA , ALPHA , DTOT 
DIMENSION SADP0T(2) 

INTEGER REGION, FORM 
REAL *8 LEMDA 

TYPE *, ’INPUT DTOT, ALPHA, S,ZL,ZT,NZ,OMEGAL,OMEGAT 
1 , NO, HDISL, HDIST, NH' 

READ *, DTOT, ALPHA, S,ZL,ZT,NZ,OMEGAL,OMEGAT, NO, HDISL, 
1 HDIST, NH 

SPEED=340. 

PI=4.*DATAN(1 .DO) 

R1 =300. 

IP(OMEGAL.EQ.OMEGAT) THEN 
STEP0=.001 

ELSE 

STEP0=( OMEGAT-OMEGAL ) /NO 
I 

DO 14 OMBGA=OMEGAL, OMEGAT , STEPO 
LEMDA=OMEGA/ (SPEED* ALPHA ) 

BRW=OMEGA/SPEED 
BRO=BRW*DSQRT (1 ./(I . +DTOT ) ) 

BRS=BRW*D3QRT ( ( 1 . +ALPHA*S ) / ( 1 . +ALPHA*S+DTOT ) ) 
ra=0MEGA/(2.*PI) 

RATI0=PR/R1 

R=1 . +9 • 08*RATI0** ( - . 75 ) 

X=-1 1 .9*RATI0**(-.73) 

Z1 =DCMPLX(R,X) 

ID=1 
ID1 =0 
ID2=0 

IF ( OMEGA. GT. 6000) THEN 
ID=3 
GOTO 1 

ENDIF 

CALL P0LE(P5,S,Z,HDIS,Z1 ,PC,ID) 

1 IF (HDIST. EQ. HDISL) THEN 
STEPH=.001 

STEPH= (HDIST-HDISL ) /NH 

ENDIF 

IF(ZT.EQ.ZL) THEN 
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STEPZ=.001 

ELSE 

STEPZ=-(ZT-ZL)/NZ 

ENDIF 

NC=0 

DO 16 Z=ZT,ZL,STEPZ 

DO 16 HDIS=HDISL, HDIST , STEPH 

CHE=BRW*HDIS 

P1=.0 

P2=.0 

P3=.0 

P4=.0 

P5=.0 

P=.0 

BRZ=BRW*DSQRT ( ( 1 . +ALPHA*Z ) / ( 1 . -(-ALPHA *Z+DTOT ) ) 
CALL SAD (S,Z,HDIS, REGION, SADPOT,ID2) 
IF(ID2.BQ.4) GOTO 2 
IF(REGION.NE.7) THEN 
BETA=SADP0T(2) 

ELSE 

BETA=SADPOT ( 1 )+SABPOT ( 2 )*DCMPLX( .0,1.) 

ENDIF 

CALL FINFORM( BETA, Z,FC,HDIS,S,F2, FORM) 

2 GOTO (10,20,10,20,30,40,50) REGION 

10 CALL INTD(SADP0T(1 ),P1,S,Z,HDIS,Z1) 

BETA=SADPGT(2) 

CALL INTR(BETA,P2,S,Z,HDIS,Z1 ,1,1 ) 

GOTO 60 

20 CALL INTD(SADP0T(1 ),P1 ,S,Z,HDIS,Z1 ) 
BETA=SADP0T(2) 

CALL INTR(EETA,P2,S,Z,HDIS,Z1 ,2,1 ) 

GOTO 60 

30 BETA=SADPOT ( 1 ) 

CALL INTR(BETA,P2,S,Z,HDIS,Z1 ,1,1 ) 
BETA=SADP0T(2) 

CALL IN7R(BETA,P1 ,S,Z,HDIS,Z1 ,2,2) 

GOTO 60 

40 IF(NC.EQ.1 ) GOTO 41 
BETA=SADPOT ( 1 ) 

CALL INTR(BETA,P2,S,Z,HDIS,Z1 ,2,1 ) 
BETA=SADP0T(2) 

CALL INTR(BETA,P1 ,S,Z,HDIS,Z' ,2,2) 

41 CALL INTC(S,Z,HDIS,Z1,P3) 

IF(NC.EQ.O) THEN 

IF(CDABS(P1 +P2).GE.CDABS(P3) ) THEN 
NC=1 
P1=.0 
P2=.0 

ELSE 

P3=.0 

ENDIF 

ENDIF 


i 1 


GOTO 60 

50 CALL INTC(S.Z,HDIS,Z1,P1) 

60 IF(IDI.EQ.I) THEN 

CALL C0N(P4,S,Z,HDIS,Z1 ) 

ELSE 

IF(ABS(FC) .LE.ABS(CHE) ) THEN 
CALL C0N(P4,S,Z,HDIS,Z1 ) 

ID1 =1 
PC -CHE 

ENDIF 

ENDIF 

61 IF(ID.EQ.O) GOTO 70 
IF(ID.EQ.I) ID=2 

IF ( OMEGA. GT. 6000) GOTO 70 
CALL P0LE(P5,S,Z,HDIS,Z1 ,FC,IDj 
70 P=P1+P2+P>P4+P5 

ID2=ID+ID1 
SPL=1 ./(4*PI*HDIS) 

SPL1 =20. *DL0G10(CDABS(P)/. 00002) 

SPL2=20 . *DL0G1 0 ( SPL/ . 00002 ) 

SPL3=SPL2-SPL1 

SPL4=20*DL0G1 0 ( CDABS ( PI +P2+P3 ) / . 00002 ) 

SPL5=20*DL0G1 0 ( CDABS ( PI +P2+P3+P4 ) / . 00002 ) 
S?L6=20*DL0G1 0 ( CDA3S ( PI +P2+P3+P4+P5 ) / . 00002 ) 
IF(Z.GT.S) THM 
BB=BRS 

EI^E 

BB=BRZ 

ENDIF 

CALL FIN(Z,S,BR0,BB,RC1 ) 

RC2=DCMPLX(RC1 ,.0D0) 

RC=CDAB3 (-G32PI ( Z , RC2)-G32PI ( S , RC2 ) ) 

IF(HDIS.GT.RC) THEN 

SPL21 =S PL2 -20 . '^DLOG 1 0 ( HD I S /RC ) / LOG 10(2. ) 
IF((SPL21-SPL2).LE.-30.) SPL21 =SPL2-30. 
SPL2=SPL21 

ENDIF 

IF(OMEGAT.NE.OMEGAL) THEN 

WRITE ( 20, *)FR,5PL1 ,SPL2,SPL3 
WRITE ( 24 , * )FR , SPL2 , SPL4 , SPL5 , SPL6 

ENDIF 

IF(ZL.NE.ZT) THEN 

WRITE(20,*)Z,SPL1 ,SPL2,SPL3 
WRITF ( 24 , * ) Z , SPL2 , SPL4 , SPL5 , SPL6 

ENDIF 

IF(HDISL.NE.HDIST) THEN 

WRITE ( 20 ,* )DL0G1 0 (HDIS ), SPL1 , SPL2 , SPL3 
WRITE( 24 , * )DL0G1 0 (HDIS ) , SPL2 , SPL4 , SPL5 , SPL6 

ENDIF 

1 6 CONTINUE 

14 CONTINUE 

STOP 


d*- 
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SUBROUTINE INTR(SADP,P,S,Z,HDIS,Z1 ,K,N) 


* THIS SUBROUTINE USES SADDLE POINT METHOD TO 

* CALCULATE THE INTEGRAL OF THE REFLECTED/REFRACTED 


* TERM. 


* 

* 

* 

* 

* 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*16 SADP, P,Q, BETA, AI,F, FI ,F2,Z1 ,E,H2,GZ 
COMPLEX *1 6 G32PI2 , C , GZZ , GZS , H2Z , H2S , HI 
COMMON /CONSTANT /S PEED , OMEGA , ALPHA , DTOT 
REAL *S LEMDA 
INTEGER FORM 
P=DCMPLX(.0, .0) 

LEMDA=OMEG A/ ( SPEED* ALPHA ) 

AI=DCMPLX(.0, 1 . ) 

PI=4.*DATAN(1 .DO) 

CALL FINFORM ( SADP, Z,FF,HDIS,S,FFF, FORM) 

C=1. 

N1 =1 
N2=1 
N3=1 

IF^FORM.EQ.3) THEN 
N3=2 

C=CDEXP(-AI*Pl/3. ) 

ENDIF 

IF(F0RM.EQ.4) THEN 
N1 =2 
N3=2 

C=CDEXP(-2.*AI*Pl/3.) 

ENDIF 

IF(F0RM.EQ.5) THEN 
N1 =2 

C=CDEXP(-AI*PI/3. ) 

ENDIF 

CALL FALL(Z,SADP,F,F1 ,F2,Z1 ,E,H2Z,GZZ,K,H1 ,N3) 

CALL FALL(S,SADP,F,F1 ,F2,Z1 ,E,H2S,GZS,K,H1 ,N1 ) 

CALL FALL(.0,SADP,F,F1 ,F2,Z1 ,E,H2,GZ,K,H1 ,N2) 
IF(K.EQ.I) THEJ 

FI =2 . *G32PI2 ( . 0, SADP )-G32PI2 ( Z , SADP )-G32PI2 (S , SADP ) 
ELSE 

FI =-G32PI2 ( Z , SADP )-G32PI2 ( S , SADP ) 

ENDIF 

F3=DSQRT ( 2 . *PI /CDABS ( FI ) ) 

IF(N.EQ.I) THEN 
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Q=PI*AI/4. 

Q=3.*PI*AI/4. 

ENDIF 

IF(DREAL(GZZ*GZS).LT..O.AND.DIMAG(GZZ*GZS).GT..O) THEN 
SS=-1. 

ELSE 

SS=1 . 

ENDIF 

P=H2S*H2Z*SS*E*F3*CDEXP (Q )*CDSQRT (2 . *SAD?/ ( PI*HDIS ) )*C 

1 *CDEXP(-AI*SADP*HDIS+AI*PI/4. )/(24-*AI 

2 *CDSQRT ( GZS*GZZ ) *F*LHIDA** ( 2 . /3 • ) ) 

RETURN 

HID 

MMWMMWMXMMMMMXMMMhUM X WMW M MMMMM U ■< M W » R MM M MUMWMMMMMMXWMyMUM M 

SUBROUTINE INTD(SADP,P,S,Z,HDIS,Z1 ) 


* * 

* THIS SUBROUTINE USES SADDLE POINT METHOD TO * 

* CALCULATE THE INTEGRAL OF THE DIRECT TERM * 

* * 


• WMWWMWWMWMWMW 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

C0MPLEX*16 P,Q,BETA,AI,F,F1 ,F2,Z1 ,E,H2,GZZ 
C0MPLEX*1 6 G32PI2,GZS,H2Z,H2S,H1S,H1Z 
COMMON /CONSTANT/SPEED, OMEGA, ALPHA, DTOT 
REAL *8 LEMDA 
P=DCMPLX( .0, .0) 

LHTOA=0MEGA/ ( SPEED*ALPHA ) 

AI=DCMPLX( .0, 1 . ) 

PI=4.*DATAN(1 .DO) 

BETA=SADP 

K=1 

CALL FALL(Z,BETA,*\F1 ,F2,Z1 ,E,H2Z,GZZ,K,H1 Z, 1 ) 

CALL FALL(S,BETA, ,\F1 ,F2,Z1 ,E,H2S,GZS,K,H1S,1 ) 
IF(Z.GT.S) THEN 

F3=DREAL ( -G32PI2 ( Z , BETA ) +G32PI2 ( S , BETA ) ) 

H2=H1S*H2Z 

ELSE 

F3=DREAL ( G32PI 2 ( Z , BETA ) -G32PI 2 ( S , BETA ) ) 

H2=H2S*H1 Z 
ENDIF 

F3=DSQRT (2. *PI/DAB3 (F3 ) ) 

Q=PI*AI/4. 

IF(DREAL(GZZ*GZS).LT..O.AND.DIMAG(GZZ*GZS).GT..O) THEN 
SS=-1. 

ELSE 

SS=1 . 
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INDIF 

P=H2*SS »F3*CDEXP ( Q)*DSQRT ( 2 . *SADP/ ( PI *HDIS ) ) 

1 *CDEXP(-AI*SADP*HDIS+AI*PI/4. )/(24-*AI 

2 *CDSQRT ( GZS*GZZ ) *LEMDA** ( 2 . /3 • ) ) 

RETURN 

IND 


SUBROUTINE INTC(S,Z,HDIS,Z1 , P) 


* THIS SUBROUTINE USES THE NONUNIFORM ASYMPTOTIC * 

* SOLUTION OF SACHS AND SILBIGER'S TO EVALUATE THE * 

* PRESSURE LEVEL NEAR THE CAUSTIC. * 

* * 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX* 1 6 P,AIQ,BETA,AI t H*,F,F1 ,F2,Z1 ,E,H2,GZZ,GZS,H2Z, 

1 H21 ,H22,H23,H11 ,H12,H13,ARG,G32PI3,BIQ, 

2 AIP,BIP,G32PI,H2S,H1 
COMMON /BRANCH/BRO , BRS , BRZ , BRW 
COMMON /CONSTANT/SPEED , OMEGA , ALPHA , DTOT 
REAL *8 LEMDA 

IF(Z.GT.S) THEN 
B1 -BRS 

ELSE 

B1 =BRZ 

ENDIF 

CALL FIN(Z,S,BRO,B1 ,B) 

BETA=B 

P=DCMPLX(.0, .0) 

LSMDA=OMEGA/ (S?EED*ALPHA ) 

AI=DCMPLX(.0,1 . ) 

PI=4.*DATAN(1 .DO) 

R1 =CDABS ( -G32PI ( Z , BET A ) -G32PI ( S , BETA ) ) 

DR=HDIS-R1 

K=1 

CALL FALL(Z f BETA,F,F1 ,F2,Z1 ,E,H2Z,GZZ,K,H1 ,1 ) 

CALL FALL(S t BETA,F,F1 ,F2,Z1 ,E,H2S,GZS,K,H1 ,1 ) 
Q43=DREAL(G32PI3 ( Z , BETA )-KJ32PI3 ( S , BETA ) ) 
ARG=(2./Q43)**(1./3.)*DR 
CALL CGBAIR(ARG,AIQ,BIQ,AIP,BIP) 

IF(DREAL(GZZ*GZS).LT..O.aND.DIMAG(GZZ*GZS).GT..O) then 
SS=-1. 

ELSE 

SS=1 . 

ENDIF 

P=H2Z*H2S*SS*DSQRT ( 2 . *B/ ( PI *HDI S ) ) 

1 *CDEXP(-AI*B*HDIS+AI*PI/4. )*ARG*2.*PI*AIQ/(24.*AI 


) 
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2 *CDSQRT (GZS*GZZ ) *LEMDA** ( 2 . /3 • ) *DR ) 

RETURN 
END 


SUBROUTINE P0LE(P5,S,Z,HDIS,Z1 , CHE, ID) 


* THIS IS A SUBROUTINE TO EVALUATE THE CONTRIBU- * 

* TION OF THE POLE. * 

* BETA IS THE LOCATION OF THE POLE IN THE COMPLEX * 

* BETA SPACE. * 

* * 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMPLEX* 1 6 AI,BETA,Z1 ,F,F1 ,F2,E,H2Z,H2S,GZZ,GZS,F1 1 ,E1 , 
1 P5,G32,H1 

COMMON /BRANCH/BRO , BR3 , BRZ , BRW 
COMMON /CONSTANT /SPEED , OMEGA , ALPHA , DTOT 
INTEGER REGION, FORM 
REAL *8 LEMDA 
jjEHDA=OMEGA/ ( 3PEED*ALPHA ) 

PI=4.*DATAN(1 .DO) 

AI=DCMPLX( .0,1.) 

N=1 
K=1 
P5=.0 

GOTO (1,20,30) ID 
1 BETA=OMEGA/SPEED*CDSQRT ( 1 .-(1 ./Z1 )**2. ) 

10 CALL FALL(0. , BETA, F, FI ,F2,Z1 ,E,K2Z,GZZ,N,H1 ,1 ) 

IF(CDABS (F).LT.. 00000000001 ) RETURN 
BETA=BETA-2 . *F/ ( FI +CDSQRT ( FI **2 . -2 . *F*F2 ) ) 

K=K+1 

IF(K.GE.500) RETURN 
GOTO 10 

20 CALL FINFORM(BETA,Z,FI’,HDIS,S,FFF,FORM) 

IF(F0RM.NE.1 ) THEN 

ID=0 

RETURN 

ENDIF 

IF(DREAL(BETA) .GT.BRW) THEN 
IF(DIMAG(BETA).GT..O.AND.ABS(FF).LE.ABS(CHE)) GOTO 30 
IF(DL4AG(BETA).LE..0.AND.ABS(FF).GE.ABS(CHE)) GOTO 30 
ELSE 

IF(DREAL(BETA).LT.BRO) THEN 

IF(DIMAG(BETA).GT..O.AND.ABS(FF).GE.ABS(CHE)) GOTO 30 
IF(DIMAG(BETA).LE..O.AND.ABS(FF).LE.ABS(CHE)) GOTO 30 
ENDIF 
ENDIF 
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RETURN 

30 CALL FALL(Z,BETA,F,F1 1 ,F2,Z1 ,E1 ,H2Z,GZZ,N,H1 ,1 ) 

CALL FALL(S,BETA,F,F1 1 ,F2,Z1 ,E1 ,H2S,GZS,N,H1 , 1 ) 

ID=3 

P5=CDEXP(-AI*BETA*HDIS+AI*Pl/4. )*H2S*H2Z*E*CDSQRT (2*BETA 
1 /(PI*HDIS))/(24.*DCMPLX(.0,.1 )*L0tDA**(2./3. ) 

1 *CDSQRT(GZZ*GZS)*F1 )*(-2*PI*AI) 

RETURN 

END 


SUBROUTINE C0N(P4,S,Z,HDIS,Z1 ) 


* SUBROUTINE CON EVALUATES THE INTEGRATION ALONG * 

* THE MODIFIED INTEGRATION PATH, WHICH IS REQUIRED TO * 

* KEEP THE INTEGRATION PATH PROM RUNNING INTO THE * 

* BRANCH POINT Bv. * 



IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

REAL*8 LPMDA 

C0MPLEX*16 P4,D,AI,F,F1 ,F2,Z1 

COMMON /CONSTANT /SPEED , OMEGA , ALPHA , DTOT 

LEHDA=OMBGA/ ( SPEED*ALPHA ) 

AI=DCM?LX(.0,1 . ) 

PI=4.*DATAN(1.D0) 

A=DSQRT ( 2 . *OMEGA/ ( PI*HDIS*SPEED*DTOT ) ) 

B=(1+ALPHA*Z+DT0T) 

C=(1+ALPHA*S+DT0T) 

D=A*(B*C)**(1 ./4. V(8*ALPHA*LPMDA*PT*(-AI*HDIS-2*SPEED/(3. 

1 ^KMEGA*DSQRT(DT0T))*(B^(3/2)+C**(3/2)-2*(1+DT0T)**(3/2)))) 

F=-AI*0MEGA*HDIS/SPEED-2*LEMDA*DSQRT (DTOT ) * ( SORT (B)-t-SQRT ( C ) 

1 -2*SQRT ( 1 +DTOT ) )+AI*PI/4 
FI = 1 +AI*Z 1 / ( 4* ( 1 +DTOT ) *LEMDA ) 

F2=-AI*Z1 *DSQRT ( DTOT/ ( 1 +DTOT ) ) 

P4=D* ( ( FI -F2 ) / ( FI +F2 ) ) *CDEXP ( F ) 

RETURN 

END 
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SUBROUTINE SAD(S, Z,HDIS, REGION, SADPOT, ID) 



* SUBROUTINE SAD IS USED TO FIND THE SADDLE POINT * 

* OF THE INTEGRATION. * 


c 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
DIMENSION SAD POT (2) 

COMMON /CONST ANT /SPEED , OMEGA , ALPHA , DTOT 
COMMON /BRANCH/ERO , BRS , BRZ , BRV 
INTEGER REGION, FORM 
CALL FINREG(S,Z, HDIS, REGION) 

GOTO (10,20,10,20,50,40,50) REGION 
10 CALL FINDRL(S,Z,HDIS,SADP0T(1 ),1 ) 

CALL FINDRL(S,Z,HDIS,SADP0T(2),2) 

RETURN 

20 CALL FINDRL(S,Z,HDIS,SADP0T(1 ),1 ) 

CALL FINRFC(S,Z,HDIS, REGION, SADPOT) 
RETURN 

30 CALL FINDRL(S,Z, HDIS, SADPOT (1 ), 2) 

CALL FINRFC(S,Z, HDIS, REGION, SADPOT) 
RETURN 

40 CALL FINRFC ( S , Z , HDIS , REGION, SADPOT ) 
RETURN 

50 IF(ID.EQ.4) RETURN 
IF(Z.GT.S) THEN 

CALL FIN(Z,S,BRO,BRS,B) 

ELSE 

CALL FIN(Z,S,BRO,BRZ,B) 

ENDIF 

CALL FINCPX(Z,S, HDIS, B, SADPOT) 

60 RETURN 
END 



SUBROUTINE FINDRL(S,Z,HDIS,SAD,I) 

C 

C UMUM U MMWWMWMMMUVM MVM MMWM yMU MMMMWMMMMMM WM M 

c * * 

C * SUBROUTINE FINDRL IS USED TO FIND THE * 

C * SADDLE POINT FOR THE DIRECT AND REFLECTED * 

C * WAVES. * 

C * * 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMPLEX* 16 G1 C , G2C , GO , G32PI 

COMMON /CONSTANT /SPEED , OMEGA , ALPHA , DTOT 

G1 =.0 

IF'J.EQ.I) THEN 
IF(Z.GT.S) THEN 
SS=S 


ELSE 


SS=Z 
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ENDIF 

FT.SF 

ss=o. 

G2=OMEGA*DSQRT ( ( 1 . +ALPHA*SS ) / ( 1 . +ALPHA*SS+DTOT ) ) /SPEED 
G1C=DCMPLX(G1 ,.0D0) 

G2C=DCMPLX(G2,.ODO) 

IF(I.EQ.I) THEN 
IF(Z.GT.S) THEN 

FI =-HD I S+DREAL ( G32 PI ( S , G 1 C )-G32PI (Z,G1 C) ) 

F2=-HDT j+DREAL ( G32PI ( S , G2C ) -G32PI ( Z , G2C ) ) 

ELSE 

FI =-HDIS-DREAL(G32PI (3 , G1 C )4G32PI (Z , G1 C ) ) 

F2=-HDIS-DREAL (G32PI (S,G2C)+G32PI(Z t G2C ) ) 

ENDIF 

ELSE 

FI =-HDIS+2 . *DREAL(G32PI ( . 0, G1 C )-G32PI ( Z , G1 C )-G32PI (S , G1 C ) ) 
F2=-HDIS+2 . *DREAL ( G32PI ( . 0 , G2C ) -432PI ( Z , G2C ) -G32PI ( S , G2C ) ) 

endif 

IF(DABS(F1 ).LT. .ID-8) GOTO 11 
IF(DAB3(F2).LT. .ID-8) GOTO 21 
IF(F1 *F2.GT. .0) GOTO 41 
K=0 

10 K=K+1 
IF(K.GT.500) GOTO 31 
G=(G1 +G2)/2. 

GC=DCMPLX(G,O.DO) 

IF(I.EQ. 1 ) THEN 

IF(Z.GT.S) THEN 

F=-HDIS+DREAL(G32PI (S, GO )^}32PI (Z , GO ) ) 

ELSE 

F=-HDIS-DREAL ( G32PI ( S , GO ) -KJ32PI ( Z , GO ) ) 

ENDIF 

ELSE 

F=-HDI S+2 . *DREAL ( G32PI ( . 0, GO )-G32PI ( Z , GO ) -G32PI ( S , GO ) ) 
ENDIF 

IF(DA33(F).LT.. ID-10) GOTO 31 
IF(F1 *F.GT. .0) THEN 
G1=G 
GOTO 10 

ELSE 

G2=G 
GOTO 10 

ENDIF 

11 SAD=G1 
F=F1 
RETURN 

21 SAD=G2 

F=F2 
RETURN 
31 SAD=G 

41 RETURN 
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SUBROUTINE FINRFC(S,Z,HDIS, REGION, SADRFC) 

n R R R W W R R R^TR R^^R R R R RR R R^r R 


* * 

* SUBROUTINE FINRPC IS USED TO FIND THE * 

* SADDLE POINTS OF THE REFRACTED TIHM. * 

* * 


yyy yyyyyyyyyyyyyHy mtMyyMyyttyyyyyyyMMyyyyyyMyMMyyy y 
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IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

INTEGER REGION 

C0MPLEX*1 6 G1C,G2C.GC,GSC,GLC,G32PI 
DIMENSION SADRFC (2) ,RFC(2) 

COMMON /CONSTANT/SPEED, OMEGA, ALPHA, MOT 
G1 =OMEGA*DSQRT ( 1 . / ( 1 . +DTOT ) ) /SPEED 
IF(Z.GT.S) THEN 
3S=S 

ELSE 

SS=Z 

ENDIF 

G2=0MEGA*DSQRT ( ( 1 . +ALPHA*SS ) / ( 1 . +ALPHA*SS+DTOT ) ) /SPEED 
G1C=DCMPLX(G1 ,6. DO) 

G2C=DCMPLX(G2,O.DO) 

FI =-HDIS-DREAL(G32PI ( Z , G1 C )-G32PI ( S , G1 C ) ) 
F2=-HDIS-DREAL ( G32PI ( Z , G2C ) -G32PI ( S , G2C ) ) 

IF (REGION. EQ. 6) GOTO 71 
IF(DABS(F1).LT..1D-8) GOTO 11 
IF(DA3S(F2).LT. .ID-8) GOTO 21 
IF(F1*F2.GT..O) GOTO 41 
K=0 

10 A=K+1 
IF(K.GT.500) GOTO 51 
G=(G1+G2)/2. 

GC=DCMPLX(G,.0D0) 

F=-HDIS-DREAL ( G32PI ( Z , GC ) -G32PI ( S , GC ) ) 

IF(DABS(F).LT.. ID-10) GOTO 31 
IF(F1 *F.GT. .0) THEN 
G1 =G 
GOTO 10 

ELSE 

G2=G 
GOTO 10 

ENDIF 

11 SADRFC (2)=G1 
F=F1 
RETURN 

21 SADRFC (2 )=G2 

F=F2 


a 
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RETURN 

31 SADRPC(2)=G 

41 RETURN 

51 SADRFC(2)=G 

RETURN 

71 IP(F1 *F2.LT. .0) GOTO 81 
CALL FIN(Z f S,G1 ,G2,G) 

GC=DCMPLX(G, .ODO) 

F=-HDIS-DREAL (G32PI ( Z , GO ) -G32PI ( S , GC ) ) 

IF(DABS(F) .LT. .ID-10) GOTO 72 
IF(F1 *F.GT. .0) GOTO 81 
IF(G2.EQ.G.0R.G1 .EQ.G) GOTO 72 
GG=G 
K=0 

20 K=K+1 

IF(K.GT.500) GOTO 73 
GS=(GG-fG1 )/2. 

GSC =DCMPLX ( GS , . ODO ) 

F3=-HDIS-DREAL ( G32PI ( Z , GSC )-G32PI ( S , GSC ) ) 
IF(DAB3(FS).LT.. ID-10) GOTO 73 
IF(FS*F1 .GT..O) THEN 
G1 =GS 

ELSE 

GG=GS 

ENDIF 
GOTO 20 

73 SADRFC(1 )=GS 
K=0 

30 K=K+1 

IF(K.GT.500) GOTO 74 
GL=(G-+G2)/2. 

GLC =DCMPLX ( GL , . ODO ) 

FL=-HDI S-DREAL ( G32PI ( Z , GLC ) -G32PI ( S , GLC ) ) 
IF(DABS(FL).LT.. ID-10) GOTO 74 
IF(FL*F2.GT. .0) THEN 
G2=GL 

ELSE 

G=GL 

ENDIF 
GOTO 30 

74 SADRFC(2)=GL 
RETURN 

72 SADRFC(1 )=G 
SADRFC ( 2 ) =G 

81 RETURN 

END 
C 
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c 

SUBROUTINE FINCPX(Z,S,HDIS,G, SADPOT) 

C 
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c * 

C * SUBROUTINE FINCPX IS USED TO FIND THE 

C * SADDLE POINTS OF A POINT INSIDE THE SHADOW 

C * REGION. 

C * 

C 

IMPLICIT REAL*8 (A-H.O-Z) 

COMMON /CONSTANT /SPEED, OMEGA, ALPHA, DTOT 
COMMON /ERANCH/BRO , BRS , BRZ , BRW 
DIMENSION SAD POT (2) 

COMPLEX* 16 G32PI , G32PI 2 , G32PI 3 , B , F , FP , FP2 
COMPLEX* 1 6 GG,G1,G2,GC 
GC=DCMPLX(G, . ODO) 

FF=-HDIS-DREAL ( G32PI ( Z , GC ) -G32PI ( S , GC ) ) 

IF(DAB3(FF).LT. .ID-6) GOTO 710 

K=0 

GG=G*DCMPLX( 1 . , . 0)+. 1 D-4*DCMPLX( .0,1.) 
FG=-HDIS-DREAL( G32PI ( Z , GG ) -KJ32PI (S,GG) ) 

11 G1=GG+.01*DCMPLX(.0,1.) 

FI =-HDIS-DREAL(G32PI(Z,G1 )-KJ32PI (S,G1 ) ) 
IF(DAB3(F1).LT..1D-6) GOTO 720 
IF(K.GT.500) GOTO 900 
IF((FF-FG)*(FF-F1 ).GT. .0) THEN 
GG=G1 
K=K+1 
GOTO 11 

. ELSE 
ENDIF 

20 G2=(GG-Kz1 )/2. 

F?=-HDIS-DREAL( G32PI ( Z , G2 )-*G32PI ( S , G2 ) ) 
IF(DABS(F2) .LT. . 1D-6.0R.K.GT.200) GOTO 730 
IF( (FF-FG)*(FF-F2).GT. .0) THEN 
GG=G2 
K=K+1 

ELSE 

G1 =G2 
K=K+1 

ENDIF 
GOTO 20 
720 GG=G1 

GOTO 900 
730 GG=G2 

GOTO 900 

710 SADP0T(1 )=G 
SADP0T(2)=.O 
RETURN 
900 B=GG 
K=0 

10 F=-HDIS-G32PI(Z,B)-G32PI(S,B) 

FP=-G32PI2 ( Z , B ) -G32PI2 (S,B) 
FP2=-G32PI3(Z,B)-G32PI3(S,B) 
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IP(Cr xj(F).LT. . ID-6) GOTO 100 
B=B-2 . *F/ ( FP+CDSQRT ( FP**2 . -2 . *F*FP2 ) ) 

K=K+1 

IP(K.GT.500) GOTO 200 
GOTO 10 

100 SADP0T(1 )=DREAL(B) 

SADP0T(2)=DIMAG(B) 

IF(DREAL(B).LE.BRS.AND.DREAL(B).LE.BRZ) 

1 SADP0T(2)=ABS(SADP0T(2)) 

200 RETURN 
END 
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SUBROUTINE FIN(Z,S,G1 ,G2,G) 



* SUBROUTINE PIN IS USED TO FIND THE * 

* BETA WHICH CORRESPONDS TO THE CAUSTIC AT * 

* RECEIVER HEIGHT Z. * 

* * 
WMWWMWMWWWWWWWWWWMMMWM W MWWMMWWWWWWMWW l MWWWWWWMWWWW 
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IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX* 1 6 G1 C , G2C , GO , G32PI2 
COMMON /CONSTANT /SPEED, OMEGA, ALPHA, DT0T 
C=.1D-6 

G1 C=DCMPLX(G1 ,0.D0) 

FI =DREAL ( G32PI2 ( Z , G1 C ) +032PI2 ( S , G1 C ) ) 

11 G2=G2-C 

G2C=DCMPLX(G2,0.D0) 

F2=DREAL ( G32PI2 ( Z , G2C ) -KJ32PI2 ( S , G2C ) ) 

IF'(DABS(F1 ).LT. .ID-8) GOTO 100 
IF(DABS(F2).LT. .ID-8) GOTO 200 
IF(F1*F2.GT..O) GOTO 400 
K=0 
GL=G1 
GR=G2 
G2=G2+C 
10 K=K+1 

IF(K.GT.200) RETURN 

G=(GL+GR)/2 

G0-DCMPLX(G,0.D0) 

F=DREAL (G32PI2 ( Z , GC )-KJ32PI2 ( S , GC ) ) 
IF(DABS(F).LT..1D-8) RETURN 
IF(P*F1.GT..0) THEN 
GL=G 

ELSE 

GR=G 

ENDIF 
GOTO 10 
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100 G=G1 
RETURN 
200 G=G2 

G2=G2+C 
RETURN 
400 C=C*.1 

GOTO 11 
END 
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SUBROUTINE FINREG; 3, Z.HDIS, REGION) 



* THIS IS A SUBROUTINE TO FIND THE REGION IN THE * 

* PHYSICAL SPACE WHERE A SPECIFIC POINT IS LOCATED. * 



IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

INTEGER REGION 

COMMON /CONSTANT /SPEED, OMEGA, ALPHA, BT0T 
THE FOLLOWING STATMENT IS USED BECAUSE OF THE 
PROBLEM ARISED FROM THE DOUBLE PRECISION. 
IF(ABS(Z-S).LT.O. IE-10) GOTO 5 
IF(Z.LE.S) THEN 
IF(Z.LT.1E-5)GOT0 8 
5 ANGLE2=DAC0S ( DSQRT ( ( 1 . +ALPHA*S+DTOT ) / ( ( 1 . +DT0T ) 

1 *( 1 .+ALPHA*S) ) ) ) 

B0UND2=FNRAY ( Z , ANGLE2 , S )+FNRAY(S, ANGLE2 , S ) 

8 IF((Z.GT.2.99999)*AND. (Z.LT.3*00001 ))THEN 

BOUND3=0.0 
GOTO 10 

ENDIF 

ANGLE3=DAC0S ( DSQRT ( ( 1 . +ALPHA*S+DTOT ) ■ * ( 1 . +ALPHA* 

1 Z)/((1.+ALPHA*Z+DT0T)*(1.+ALPHA*S)))) 

B0UND3=FNRAY ( S , ANGLE3 , S ) 

IF(Z.LT. 1 E-5 ) B0UND2 =B0UND3 

10 IF((HDIS.GT.B0UND2) .AND. (HDIS.GT.B0UND3) ) THEN 

CALL EVALB4(S,Z,ANGLE4,B0UND4) 
IF(HDIS.LE.B0UND4) THEN 
REGI0N=6 
ELSE 

REGI0N=7 

ENDIF 

ELSE 

IF( (HDIS.GT.B0UND2).AND. (HDIS.LT.B0UND3)) THEN 
REGION =4 

EL3EIF( (HDIS.GT.B0UND3)«AND. (HDIS.LT.B0UND2) ) THEN 
REGI0N=5 


ELSE 
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► REGION =3 

ENDIF 

ENDIF 

ZER0=0.0 

B0UND1 rFNRAY(Z,ZERO,S) 

> ANGLE2=DAC0S ( DSQRT ( ( 1 . + ALPHA *S+DTOT ) / ( ( 1 . +DT0T ) 

1 *(1 .+ALPHA*S)')) 

B0UND2=FNRAY ( Z , ANGLE2 , S)+FNRAY(S, ANGLE2 . S ) 

IF( (HDIS. LT.B0UND2). AND. (HDIS.LT.D0UND1 )) THEN 
REGI0N=1 

ELSEIF((HDIS.LT.B0UND1 ).AND. (HDIS.GT.B0UNI)2)) THEN 
REGION =2 

ELSEIF( ( HDIS. LT.B0UND2). AND. (HDIS. GT. BOUND 1 )) THEN 
REGI0N=5 

ELSE 

CALL EVALB4 ( S , Z , ANGLE4 , B0UND4 ) 

IF (HDIS . LT . B0UND4 ) THEN 
REGION =6 
GOTO 20 
ENDIF 
REG10N=7 
ENDIF 
ENDIF 
20 RETURN 
END 


SUBROUTINE EVALB4(S,Z,ANGLE4,B0UND4) 

w M w w y yM y MMMMM y MyMyyy y MyyMMMy w wyywy yww ywuyyMyyyyywu 

* * 


* SUBROUTINE EVALB4 IS USED TO EVALUATE * 

* THE SHADOW BOUNDARY. * 



IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
C0MPLEX*1 6 GUI ,GU2,G32PI2 
COMMON /CONSTANT /SPEED , OMEGA , ALPHA , DTOT 
GUESS2= (OMEGA/SPEED )*DSQRT(1 ./(I . +DT0T ) ) 
IF(Z.GT.S) THEN 

GUESS1 = (OMEGA/ SPEED )*DSQRT((1 .+ALPHA*S)/ 
1 ( 1 . +ALPHA*S+DTOT ) )-0. 00001 

GUESS 1 = ( OMBGA/SPEED ) *DSQRT ( ( 1 . -rALPHA*Z ) / 
1 ( 1 . +ALPHA*Z+DTOT ) )-0. 00001 

ENDIF 
IC0UNT=0 
10 GUI =DCMPLX( GUESS 1 ,0.D0) 

Y=DREAL(G32PI2(Z,GU1 )+G32PI2(S,GU1 ) ) 
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IF(Y.LT.O.O) THEN 
GOTO 20 
ELSE 
TB*P=Y 

IF(Z.LT.25.0) THEN 

GUESS 1 =GUESS1 -.001 

ELSE 

GUESS 1 =GUESS1 0001 

END IF 
GOTO 10 
END IF 

20 IF(Z.LT.25.0) THEN 

GUESS 1 =GUESS1 +.001 
ELSE 

GUESS 1 =GUESS1 + .0001 
END IF 
FF=TEMP 

50 GU2=DCMPLX(GUESS2,O.DO) 

HDIS=DREA1 ( G32PI2 ( Z , GU2 )+032PI2 ( S , GU2 ) ) 

GG=HDIS 

DELTA=GUESS2-(GUESS2-GUESS1 )/(GG-FF)*GG 
IF (ABS(DELTA-GUESS2).LT.. 000001 ) THEN 
GOTO 40 
ELSE 

IC0UNT=IC0UNT+1 
IF(ICOUNT.GT.IOOO) THEN 
GOTO 50 

ELSE 

GUESS2=DELTA 
GOTO 30 

END IF 
END IF 

40 ABETA=DELTA 

ANGLE4=DAC0S ( (ABETA*3PEED ) / ( OMEGA*DSQRT ( ( 1 . + 

1 ALPHA # S ) / ( 1 . +ALPHA*S+DTOT ) ) ) ) 

30UND4=FNRAY ( Z , ANGLE4 , S )+FNRAY ( S , ANGLE4 , S ) 

50 RETURN 

END 
C 

0 »» »»»»»»»»«-»■» »» » X H l> I K I I It » » »»»»»»»» 

c 

SUBROUTINE FALL(Z,BETA,F,F1 ,F2,Z1 ,E,H2,GZ,K,H1 ,N) 

C 

C * * 

C * THIS SUBROUTINE EVALUATES THE REFLECTED/REFRA- * 

C * CTED COEFFICIENTS AND THE FIRST AND SECOND DERIVA- * 

C * TIVES OF THE DENOMINATOR OF THE COEFFICIENT. * 

C * * 

C 

IMPLICIT REAL*8(A-H,0-Z) 
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COMPLEX *16 BETA, F, FI ,F2,A,GZ,GZZ,GZ1 ,GZ2,GZZ1 t GZZ2 

1 Hi, Ell ,EI2,A1 ,A2,AI,Z1 ,B,B1 ,B2,E,H1 ,H1 1 ,H12,H2,H21 

2 ,II13,C,AIQ,H22,H23 

COMMON / CONSTANT/SPEED , OMEGA , ALPHA , DTOT 

CALL GZALL1 (Z,EETA,GZ,GZ1 ,GZ2,GZZ,GZZ1 ,GZZ2,EN,Hi1 ,Hi2) 

PI=4.*DATAN(1 .DO) 

AI=DCMPLX(0. ,1 . ) 

IF(N.EQ.2) Ei=EN*CDEXP(2.*AI*PI/3.) 

CALL HALL(EN,H2,H21 ,H22,H23,H1 ,H1 1 ,H12,H13,AIQ) 

A=1 . -AI*SPEED*Z1 *GZZ/ ( 2 . *OKEGA*GZ ) 

A1 = ( A— 1 . )*GZZ1 /GZZ+(1 .-A)*GZ1 /GZ 

A2=(A-1 . )*GZZ2/GZZ+(1 .-A)*GZ2/GZ~2.*A1 *GZ1 /GZ 

B=AI*3PEED*Z 1 * ( 3 . *OMBGA/ ( 2 . *SPEED*ALPHA ) ) ** ( 2 . /3 . ) *GZ /OMEGA 

B1 =B*GZ1 /GZ 

B2=B*GZ2/GZ 

IP(K.EQ.2) GOTO 10 

P=A*H2+B*H21 

E=-(A*H1+B*H1 1 ) 

FI =A*H21 *EN1 +H2*A1 +B*H22*EN1 +B1 *H21 
F2=(A*H22+B*H23)*EN1 **2+(A*H21 +B*H22)*EN2+ 

1 (H21 *A1 +H22*B1 ) >2 . *EN1 +H2*A2+H21 *B2 

RETURN 

10 PI=4.*DATAN(1 .DO) 

C=CDEXP(PI*Al/3- ) 

F=A*(H2+C*H1 )+B*(H21 +C*H1 1 ) 

E=CDEXP(-PI*Al/3. )*(A*H2+B*H21 ) 

F1=A*(H2l4C*H11 )*EN1+(H2+C*H1 )*A1+B*(H224€*H12)*EN1+B1 
1 *(H21+C*H11 ) 

F2=(A*(H22+C*H12)+B*(H23+C*H13) )*EN1 **2+(A*(H21 +C*H1 1 ) 

1 j -r jt( ’H22+C*H1 2) )*EN2+( (H21 +C*H1 1 )*A1+(H22+C*H1 2)*B1 )*2.*EN1 

2 +(u2+C*H1 )*A2+(H21 40*H1 1 )*B2 
RETURN 
EiD 
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SUBROUTINE HALL(Z,H2,H21 ,H22,H23,H1 ,H1 1 ,H12,H13,AI) 


* HALL USES SUBROUTINE CGBAIR TO CALCULATE 1/3 ORDER * 

* HANKLE FUNCTIONS FROM AIRY FUNCTIONS. * 


IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX* 16 Z,AI,BI,AIP,BIP,K,KS,H1,H2,H11 ,H21 ,ARG,CI,H12 
COMPLEX* 16 H23,H13,H22 
CI= DCMPLX(0.D0,1.D0) 

PI= 3.141592654D0 

ARG= DCMPLX(0.D0,-PI/6.D0) 

K= (12.D0)**(1 .D0/6.D0)*CDEXP(ARG) 

KS= DCONJG(K) 


CALL CGBAIR(-Z,AI,BI,AIP,BIP) 
H1= K*(AI-CI*BI) 

H 2= KS*(AI4CI*BI) 

H1 1= -K*(AIP-CI*BIP) 

H21 = -KS* ( AIP+CI*BIP ) 

H22=-Z*H2 

HI 2=^Z*H1 

H23=-H2-Z*H21 

HI 3=— HI -Z*H1 1 

RETURN 

ETO 


SUBROUTINE CGBALR(Z,AI,BI,AIP,BIP) 



* SUBROUTINE CGBAIR EVALUATES THE AIRY FUNCTIONS. * 

* * 


IMPLICIT REAL *8 (A-H.O-Z) 

COM?LEX*1 6 Z,AI,BI,AIP,BIP,ZETA,CZETA,Z14,SUM1 ,SUM2, 

1 SUM3 , SUM4 , ZETAP, FACT 1 ,FACT2,SN,CS,FTERM,FPTERM,GTERM, 

2 GPTERM,F,FP,G,GP,Z3 
DIMENSION C (21 ) ,D(21 ) 

DATA Cl ,C2,PIRT,PI4/.3550280539DO,. 25881 94038DO, 

+ 1.772453851 DO,. 7853981 635DO/ 

DATA C/1 .DO, .069444444444444D0, 

+ .0371 33487654321 DO,. 0379930591 27800D0, 

1 . 0576491 904 12669DO,. 11 60990640255 IDO, 

+ .291 591 39923074DO, .87766696950998DO, 

2 3* 0794530301 73 1 DO, 1 2 . 341 573332345DO, 

+ 55- 62278536591 4D0, 278. 46508077759D0, 

3 1533.1 6943201 27D0, 9207. 2065997258D0, 

+ 59892. 51 3565875DO, 41 9524.8751 1653DO, 

4 31 48257. 41 78666D0, 251 9891 9-871 601 DO, 

+ 21 4288036. 96366DO, 1 929375549- 1 823D0, 

5 1 8335766937. 889D0/ 

DATA D/1. DO, -.097222222222221 DO, 

+ -.0438850308641 97D0,-.042462830789894D0, 

1 -. 0626621 63492031 DO, -. 1 241 0589602727D0, -. 308253764901 07D 0, 

2 -. 92047999241 291 DO, -3. 21 04935846485D0 , -1 2.807293080735D0, 

3 -57. 50830351 391 IDO, -287. 0332^71 0920D0,-1 576. 3573033370D0, 

4 - 9446 . 3548230953D0 , - 6 1 335 • 706663847DG , -428952 . 40040004D0 , 

5 -321 4536 . 521 4006D0, -25697908. 383909D0, -21 8293420. 8321 4D0, 

6 -1 963523788. 9909D0.-1 8643931 088. 105D0/ 
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ABSZ=ABS(Z) 

IP(ABSZ.EQ.O) GOTO 3 

IF(ABS(DIMAG(Z)).LE.1.D-12.AND.DREAL(Z).LT.O.DO) GOTO 5 
ARGZ=ATAN2 ( DIMAG ( Z ) , DREAL ( Z ) ) 

GOTO 4 

3 ARGZ=O.DO 
GOTO 4 

5 ARGZ=3* 1 41 5926535898DO 

4 CONTINUE 
IF(A3SZ.GT.6.D0) GOTO 10 

C 

C ASCENDING SERIES 

C EQS. 10.4.2,10.4-3 

C 

Z3=Z**3 
FTERM=1 .DO 
PPTERM=Z*Z/ 2 . DO 
GTERM=Z 
GPTERM=1 .DO 
GLIM=1.D-1 3*ABSZ 
P=FTERM 
FP=F?TERM 
G=GTESM 
GP=GPTIRM 
DO 1 1=1 ,100 
13=3*1 

FTERM=FTERM*Z3/((I3-1 .D0)*I3) 

FPTERM=FPTERM*Z3 / ( 13* ( 13+2 . DO ) ) 

GTERM=GTERM*Z3/(I3*(I3+1 -DO)) 

GPTERM=GPTERM*Z3/ ( ( 13-2 . DO ) *13 ) 

F=F+FTERM 

FP=FP+FPTERM 

G=G+GTERM 

GP=GP+GPTERM 

IF(ABS(GTERM) .LE.GLIM) GOTO 2 

1 CONTINUE 

2 AI=C1 *F-C2*G 
AIP=C1 *FP-C2+GP 

BI=1 . 732O5O0O8DO* ( Cl *F+C2*G ) 

BIP=1 . 732050808D0* ( C 1 *FP+C2*GP ) 

GOTO 999 
C 

C ASYMPTOTIC EXPANSION FOR /Z/ LARGE 

C 

10 SIGN=1.D0 

SUM1 =0.D0 
SUM2=0.D0 
SUM3=0.D0 
SUM4=0.D0 

IF(ABS(ARGZ).GE.1 .3D0) GOTO 20 
C 

C /ARG(Z)/ LE PI/3 
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C EQS. 10.4.59, 10.4.61, 10.4.63, IO. 4.66 

C 

ZETA=CZETA ( ABSZ , ARGZ ) 

DO 11 1=1,12 
K=I-1 

ZETAP=ZETA**K 
SUM1 =SUM1 +SIGN*C ( I ) /ZETAP 
SUM2=SUM2+SIGN*D ( I ) /ZETAP 
SUM3=SUM3+C(I)/ZETAP 
SUM4=SUM4+D ( I ) / ZETAP 
1 1 SIGN=-SIGN 

Z1 4=ABSZ** . 25DO*DCMPLX (COS ( ARGZ/4 . DO) , SIN ( ARGZ/4 . DO ) ) 
FACT1 = . 5D0*EXP (-ZETA ) / ( PIRT*Z1 4 ) 

FACT2= . 5D0*EXP ( -ZET A ) *Z 1 4/PIRT 
AI=PACT1 *SUM1 
AIP=-FACT2*SUM 2 
PACT1 =EXP(ZETA)/(PIRI*Z14) 

FACT2=EXP (ZETA )*Z1 4/PIRT 
BI=PACT1 *SUM3 
BIP=FACT2*SUM4 
GOTO 999 
C 

C /ARG(Z)/ GT PI/3 

C EQS. 10.4.60, 10.4.62, 10.4.64, 10.4.67 

C 

20 ARGZ=ATAN2 ( -DIMAG ( Z ) , -DREAL ( Z ) ) 

ZET A=CZETA ( ABSZ , ARGZ ) 

DO 21 1=1 ,10 
K2=(I-1 )*2 
J=K2+1 

ZETAP=ZETA**X2 

SUM1 =SUM1 +SIGN*C(J ) /ZETAP 

SUM2=SUM2+SIGN*C ( J+1 )/(ZETAP*ZETA) 

SUM3--SUM3+SIGN*D ( J ) /ZETAP 
SUM4=SUM4+SIGN*D(J+1 )/(ZETAP*ZETA) 

21 SIGN=— SIGN 

Z 1 4 =ABSZ** . 25DO*DCMPLX ( COS ( ARGZ/4 . DO ) , SIN ( ARGZ /4 . DO ) ) 
FACT1 =1 .D0/(PIRT*Z14) 

FACT2=Z14/riRT 
SN=SIN ( ZETA+PI4 ) 

CS=C0S(ZETA+PI4) 

AI =FACT 1 * ( SN*SUM1 -CS*SUM2 ) 

AIP=-FACT2 # ( CS*SUM>+SN*SUM4 ) 

BI=FACT1 *(CS*SUM1 +SN«SUM2) 

BIP=FACT2* ( 3N*SUM3-CS*SUM4 ) 

959 RETURN 

END 
C 


► 


FJNCTION CZET A ( ABSZ , ARGZ ) 
IMPLICIT REAL*8 (A-H,0-Z) 


oooooooo ooo 
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C0MPLEX*16 CZETA 
ARG=ARGZ*1 . 5D0 

CZETA= ( ABSZ**1 . 5D0 ) *DCMPLX ( COS ( ARG ) , SIN ( ARG ) ) 
1 *. 66666666666667DO 

RETURN 
END 


SUBROUTINE FINF0RM(BETA,Z,F1 ,HDIS,S,F2,F0RH) 


* * 

* SUBROUTINE FINFORM DECIDES THE FORM OF * 

* THE SOLUTION TO BE USED. * 

* * 


IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX* 1 6 BETA,G320,G32S,G32Z,G32,G321 ,F 
INTEGER FORM 
REAL *8 LEMDA 

COMMON /BRANCH/BRO , BRS , BRZ , BRW 
COMMON /CONSTANT /SPEED , OMEGA , ALPHA , DTOT 
PI=4.*DATAN(1 .00) 

LEMDA=OMEJGA/ (SPEED* ALPHA ) 

G320=G32 ( • 0 , BETA ) 

G32S=G32 ( S , BETA ) 

G32Z=G32 ( Z , BETA ) 

FORM=1 

IF (DIMAG (BETA) .GT. .0) THEN 

IF ( DREAL ( G320 ) . LT . . 0. AND . DIMAG ( G320 ) . GT . . 0) F0RM=2 
IF(Z.GT.S) THEN 

IF ( DREAL ( G32S ) . LT . . 0 . AND . DIMAG ( G32S ) . GT . . 0 ) F0RM=5 
IF ( DREAL ( G32Z ) .LT. . 0. AND. DIMAG (G32Z ) .GT. .0) F0RM=4 
ELSE 

IF ( DREAL’ ( G32 Z ) . LT . . 0 . AND . DIMAG ( G32Z ) .GT. .0) F0RI4=3 
IF ( DREAL (G32S ) .LT. . 0. AND. DIMAG (G32S) .GT. .0) F0RM=4 
MDIF 
EI^E 

IF ( DREAL ( BETA ).LT. BRW) F0RM=4 
IF(Z.GT.S) THEN 

IF(DREAL(BETA).LT.BRZ) F0RM=5 
IF ( DREAL ( BETA ).LT. BRS) F0RM=2 
ELSE 

IF(DREAL(HETA) .LT.BRS) FORM-3 
IF ( DREAL ( BETA ).LT. BRS) FORM-2 
ENDIF 

IF(DREAL(BETA).LT.BRO) FORM-1 
ENDIF 

GOTO (10,20,30,40,50) FORM 

F= ( -BETA*HDI3+LEMDA* ( 2 . *G321 ( . 0 , BETA , S , FORM ) 
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1 -G321 (Z,BETA,S,F0RM)-G321 (S,BETA,S,FORM))) 

GOTO 55 

20 F=(-BETA*HDIS-LEMDA*(G321 (Z, BETA, S, FORM) 

1 +G321 (S, BETA, S, FORM))) 

GOTO 55 

30 F=(-BETA*HDIS-LEMDA*(G321 (Z. BETA, S, FORM) 

1 -*G321 (S, BETA, S, FORM)) -PI/3.) 

GOTO 55 

40 F= ( -BETA*HDIS-LEMDA* ( G321 ( Z , BETA , S , FORM ) 

1 -HJ321 (S, BETA, S, FORM) )-2.*PI/3.) 

GOTO 55 

50 F= ( -BET A*hD I S-LEKDA* ( G32 1 (Z, BETA, S, FORM) 

1 +G321 (S,BETA,S,FQRM))-PI/3.) 

GOTO 55 

51 F=(-EETA*HDIS-LEMDA*(G321 (Z, BETA, S, FORM) 

1 4CDABS(G321(.0, BETA. S, FORM)) 

2 -KJ321 (S, BETA, S, FORM) )+PI/4.) 

55 FI =DREAL(F) 

F2=DIMAG(F) 

60 RETURN 
END 

C 

C 

SUBROUTINE GZALL1 (Z,BETA,GZ,GZ1 ,GZ2,GZZ,GZZ1 ,GZZ2,M,M1 ,M2) 
C 

c * * 

C * GZALL1 CALCULATES ALL THE PARTIAL DERIVATIVES * 

C * OF THE G FUNCTION WITH RESPECT TO Z. * 

C * 

C 

IMPLICIT REAL*8 ( A-H , 0-Z ) 

COMPLEX *16 BETA,GZ,GZ1 ,GZ2,GZZ,GZZ1 ,GZZ2,G,GB,GEB,SQRT1 , 

1 IN,EN1 ,IH2 

COMMON /CONSTANT/SPEED , OMEGA , ALPHA , DTOT 
CALL GALL(Z , BETA , G, GB, GBB ) 

A=1 . +ALPHA*Z+DT0T 
BRAN=OMEGA/SPEED*DSQRT ( ( A-DTOT ) /A ) 

IF(DREAL(G) .LE. .0. AND.DIMAG(G) .GE.O. ) THEN 
SI=*-1 . 

ELSE 

SI=1 . 

END IF 

GZ=SI*2 . *ALPHA*SQRT 1 ( BETA , Z ) / ( 3 . *CDSQRT (G*A) ) 

C=2 . *ALPHA**3*DT0T/ ( 9*A**2 ) 

GZZ =C/ ( GZ*G )- . 5*GZ**2/G 
B=-(2.*ALPHA*SPEED/ (3»*0MEGA) )**2 
GZ1 =B*BETA/(GZ*G)~. 5*GZ*GB/G 

GZ2=B* ( GZ*G-EETA* ( GZ 1 «G40Z*GB ) ) / ( GZ*G ) **2- . 5* ( G* ( GZ1 *GB+ 

1 GZ*GBB ) -GZ*GB**2 ) /G**2 


oooooooo ooo oooooooo ooo 
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GZZ 1 =-C* ( GZ*GB+GZ 1*G)/(GZ*G) **2- . 5* ( 2 . *G*GZ*GZ 1 -GZ**2*GB ) 

2 /G**2 

GZZ2=-C* ( GZ*G* ( GZ*GBB+2 . *GZ1 *GB+G*GZ2 )-2 . *(GZ*GB+GZ1 *G)**2) 

3 / (GZ*G ) **3-. 5* ( G* ( 2 . *G*GZ1 **2+2 . *G*GZ*GZ2-GZ**2*GBB)-2 . *GB* 

4 ( 2 . *G*GZ*GZ 1 -GZ**2*GB ) ) / G**3 
T=(3.*OMBGA/(2.*ALPHA*SPEED))**(2./3.) 

Bf=T*G 

Ei1=T*GB 

EN2=T*GBB 

RETURN 

END 


SUBROUTINE GALL(Z, BETA, G, GPI, GPI2) 


* * 

* GALL EVALUATES TEE G FUNCTION AND ITS * 

* PARTIAL DERIVATIVES WITH RESPECT TO BETA. * 

* * 


LMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX *16 G, GPI , GPI2 , BETA , G32 , G32PI , G32PI2 
COMMON / CONST ANT/SPEED , OMEGA , ALPHA , DTOT 
G=CDEXP ( 2 . /3 . *CDL0G ( G32 ( Z , BETA ) ) ) 
I?(DRSAL(G).LE..O.AND.DIMAG(G).GE..O) THEN 
S=-1. 

M£E 

S=1 . 

ENDIF 

GPI =2 . *SPEED* ALPHA *S*G32PI ( Z , BETA ) / ( 3 . *OMEGA*CDSQRT (G) ) 
GPI2=GPI*G32PI2 ( Z , BETA ) /G32PI (Z , BETA)-. 5*GPI**2/G 
RETURN 
END 
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FUNCTION G321 (Z, BETA, S, FORM) 



* THIS FUNCTION IS REQUIRED BECAUSE OF DIFFERENT * 

* BRANCH CUTS INVOLVED WITH HANKEL FUNCTIONS. * 



IMPLICIT REAL*8 (A-H.O-Z) 
COMPLEX* 1 6 BETA.G321 ,G,G32,C 
INTEGER FORM 
PI=4-*DATAN(1 .DO) 


OOOOOOOOO OOOOOOOOOOO 


85 


C=CDEXP(2.*PI*DCMPLX(.0,1 .)/3.) 

G^DEXP(2./3.*CDL0G(G32(Z,BETA))) 

G321=(CDSQRT(G))**3. 

I P ( FORM . EQ . 3 . AND . Z . LT . S . OR . FORM . BQ . 5 • AND . Z . EQ . S 
1 .0R.F0RM.EQ.4) THEN 

IF(DIMAG(BETA ) . LT . . 0)THEN 

G321 =-(CDSQRT(C*G) )**3. 

ELSE 

G321=(CDSQRT(C*G))**3. 

ENDIF 

ENDIF 

RETURN 

END 


FUNCTION G32(Z,BETA) 


* THIS FUNCTION CALCULATES THE G TO THE * 

* 3/2 POWER. * 

* * 


IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX* 1 6 BETA, PHI, SORT 1 ,SQRT2,FOO,G32,MODLOG 
COMPLEX* 1 6 SI SB 

COMMON /CONSTANT /SPEED , OMEGA , ALPHA , DTOT 

AA=1 .+ALPHA*Z+DTOT 

B=1 . - ( SPEED*BET A/ OMEGA ) **2 

SI =SQRT1 (BETA,Z) 

S=SQRT2(BETA) 

Pffl=Sl/(S*DSQRT(AA)) 

IF(DT0T.EQ. .0) THEN 
F00=.0 

ELSE 

FOO=MODLOG( (1 .+PHI)/(1 .-PHI)) 

ENDIF 

G32=DSQRT ( AA ) *31 -. 5*DT0T*F00/S 

RETURN 

END 

FUNCTION G32PI(Z,EETA) 



* THE FIRST DERIVATIVE OF G32 WITH RESPECT * 

* TO BETA IS CALCULATED. * 

* * 


oooooooo ooo 
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IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEXES BET A, PHI, SORT 1 ,SQRT2,F00,M0DL0G,G32PI,S1 ,S,A,B 

COMMON /CONSTANT/SPEED, OMEGA, ALPHA, DTOT 

AA=1 . +ALPHA*Z+DTOT 

B=1 . - ( SPEED*BET A/ OMEGA ) **2 

SUSQRT1 (BETA,Z) 

S=SQRT2(BETA) 

PHI=S1/(3*DSQRT(AA)) 

IF(DTOT.EQ. .0) THEN 
F00=.0 

ELSE 

FOO=MODLOG ( ( 1 .+PHI)/(1 .-PHI)) 

ENDIF 

A=-1 . *SPEED*BET A/ ( ALPHA*OMEGA*B ) 
G32PI=A*(DSQRT(AA)*S1+.5*DTOT*FOO/S) 

RETURN 

END 


l/WW M MW W WMMMMM X V M W W U 


FUNCTION G32PI2(Z,BETA) 


WWWMMMMMXMMMMMWMM.'KW M _ UtfWVyWWMMMMUUVMUWWUMMVMWMMWM 

* * 

* THE SECOND DERIVATIVE OF G32 WITH * 

* RESPECT TO BETA IS CALCULATED. * 

* * 


IMPLICIT REAL*8 (A-H,0--Z) 

COMPLEX* 1 6 BETA , PHI , SQRT 1 , SQRT2 , F00 , M0DL0G, PRT2 , PRT3 

C0MPLEX*1 6 G32PI2,S1 ,S,A,B,B1 ,A0,A1 ,A2,A3,PRT0,PRT1 

COMMON /CONSTANT/SPEED, OMEGA, ALPHA, DTOT 

AA=1 . +ALPHA*Z+DTOT 

B=1 . - ( SPEED*BET A/ OMEGA ) **2 

SI =SQRT1 (BETA,Z) 

S=SQRT2(BETA) 

PHI=S1 / (S*DSQRT(AA) ) 

IF(DT0T.EQ. .0) THEN 
F00=.0 


F00=M0DL0G((1 .+PHI)/(1 .-PHI)) 

ENDIF' 


1 


B1 =1 .+2.*(3PEED*BETA/0MEGA)«-*2. 

G32PI2=-SPEED/(ALPHA*OMEGA*B**2)*( (B*AA**2-B1 *DT0T*AA)/ 
(DSQRT(AA)*S1 )+. 5*DT0T*B1 *FOO/S) 

RETURN 

END 


\ 


C 

c 
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C 

FUNCTION G32PI3(Z,BETA) 

^in w nn n w w w w w inrw w n icwinnnnr 

* * 

* THE THIRD DERIVATIVE OF G32 WITH RESPECT * 

* TO BETA IS CALCULATED. * 

* * 


IMPLICIT REALMS (A-H,0-Z) 

COMPLEX* 1 6 BET A, PHI, SORT 1 ,SQRT2,FOO,MODLOG,G32PI2,PRT3 

COMPLEX* 16 G32PI3,S1 ,S,A,E,B1 ,A0,A1 ,A2,A3,PRT0,PRI1 ,PRT2 

COITION /CONST ANT/SFEED, OMEGA, ALPHA, DTOT 

AA=1 . +ALPHA*Z+DTOT 

B=1 . - ( SPEED*3ETA/ OMEGA ) **2 

SI =SQRT1 (BETA,Z) 

S=SQRT2(BETA) 

PHI=S1 / (S*DSQRT(AA) ) 

IF(DT0T.EQ. .0) THEN 
F00=.0 

ELSE 

FOO=MODLOG( ( 1 .+FHI )/( 1 .-PHI)) 

END IF 

A1 =( SPEED*BET A /OMEGA ) **2 
A0=A1 /BETA 
A2=1 .+2.*A1 
A3=5.-2.*A1 

PRT0=- ( SPEED/ ( ALPHA*OMEGA ) ) /B**2 

PRT1 =A0*AA**2*(2.*A1 *DT0T*AA+4. *DT0T**2-B*AA**2-3.*DT0T*AA) 
PRT2=(B*AA**2-DT0T*AA)*31 *DSQRT ( AA ) 

PRT3=FOO*AO*A3/B*S-2 . *A0*A2*DSQRT ( AA ) / 

1 (B*S1 ) 

PRT4=4 . *AO/B*G32PI2 ( Z , BETA ) 

G32PI3=PRT0* ( PRT1 /PRT2+0. 5*DTOT*PRT3 )+PRT4 

RETURN 

END 

m x y w y yyMyywy y yMMyyyyyyyyyyyyyyMyyy yyvM yyyyyyuyyMMyyyyyyyu u 

T* W n nn n W W n W r A a W n ^ln Wn a a W a a a A W ^I^^a n ^T^Ia^Ia a a 


FUNCTION SQRT1 (BETA,Z) 

IMPLICIT REAL*8 (A-H,0-Z) 

C0MPLEX*16 BETA,AA,SQRT1 

COMMON /CONSTANT /SPEED , OMEGA , ALPHA , DTOT 

AA=1 .-(( SPEED/OMEGA )*BETA)**2. 

BB=1 .+AI.PHA*Z+DT0T 

BRAN1 = ( OMEGA/SPEED ) *DSQRT ( ( 1 . +ALFHA*Z ) / ( 1 . -rALTHA*Z+DT0T ) ) 
IF( (DREAL(BETA) .GE. BRAN1 ) .AND. 

1 (DIMAG (BETA) .LE. 0.0)) THEN 
SQRT1 =-CDSQRT ( AA*BB-DT0T ) 

ELSE 


SORT 1 =CDSQRT ( AA*BB-DT0T ) 


ooo ooo ooo 
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HTO IP 
RETURN 
END 

FUNCTION SQRT2(BETA) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEXES BETA, AA.SQRT2 
COMMON /BRANCH/BRO , BRS , BRZ , BRW 
COMMON /CONSTANT/SPEED , OMEGA , ALPHA , DTOT 
AA=1 .-( ( SPEED/ OMEGA ) *BETA ) ** 2 . 

IP( (DREAL(BETA) .GE. ERW) .AND. 

1 (DIMAG (BETA) .LE..O)) THEN 

SQRT2=-CDSQRT ( AA ) 

P j ^ 

SQRT2=CDSQRT ( AA ) 

END IP 
RETURN 
END 


FUNCTION MODLOG(QUAN) 

IMPLICIT REAL*8 (A-H.O-Z) 

C0MPLEX*16 QUAN,M0DL0G 

IP( (DREAL(QUAN) .LE. 0.0) .AND. ( DIMAG (QUAN) 

1 .GE. 0.0)) THEN 

M0DL0G=CDL0G ( QUAN ) +DCMPLX (0.0, -2*3 -141 5927) 

FTffE 

M0DL0G=CDL0G ( QUAN ) 

END IP 
RETURN 
END 

yMUWMyMMyMMMyyMUWMMMMMMMMyUMMMMMWWMUMy VM MWUMMMUUUMyyyyUMMW 
Wn n "n W n"W n W n W n n A 

FUNCTION FNRAY(Z, ANGLE, S) 

IMPLICIT REALMS (A-H,0-Z) 

COMPLEX* 1 6 G32PI,BC 

COMMON /CONSTANT/SPEED, OMEGA, ALPHA, DTOT 
IF(ANGLE.EQ. .0) ANGLE=. 000001 

BETA=0MEGA/SPEED*DSQRT ( ( 1 . +ALPHA*S ) / ( 1 . +ALPHA*S+DT0T ) ) 

1 *DC0S( ANGLE) 

BC=DCMPLX(BETA, .0D0) 

PNRAY=-DREAL ( G32PI ( Z , BC ) ) 

RETURN 

END 


REFERENCES 


1. U. Ingard, "On the Reflection of a Sperhical Sound Wave From an 
Infinite plane," J. Acoust. Soc. Am. , Vol . 23, 1951 , pp. 329-335. 

2. R.B. Lawhead and I. Rundick, "Measurements on an Acoustic Wave Pro- 
pagated Along a Boundary," J. Acoust. Soc. Am., Vol. 23, 1951, pp. 
541-545. 

3. R.B. Lawhead and I. Rundick, "Acoustic Wave Propagation Along a 
Constant Normal Impedance Boundary," J. Acoust. Soc. Am. , Vol. 23, 
1951, pp. 546-549. 

4. A.R. Wenzel, "Propagation of Waves Along an Impedance Boundary," 
J. Acoust. Soc. Am. , Vol. 55, 1974, pp. 956-963. 

5. C.F. Chien and W.W. Soroka, "Sound Propagation Along an Impedance 
Plane," J. Sound Vib. , Vol. 43, 1975, pp. 9-20. 

6. S. Thomasson, "Reflection of Waves From a Point Source by an Impe- 
dance Boundary," J. Acoust. Soc. Am. , Vol. 59, 1976, pp. 780-785. 

7. T.F. Emblet.on, J.E. Piercy and N. Olson, "Effect of the Ground on 
Near-horizontal Sound Propagation," J. Acoust. Soc. Am., Vol. 53, 
1973, pg. 340(A). 

8. T.F. Embleton, J.E. Piercy and N. Olson, "Outdoor Sound Progaga- 
tion Over Ground of Finite Impedance," J. Acoust. Soc. Am., Vol. 
59, 1976, pp. 267-277. 

9. F.M. Wiener and D.N. Keast, "Experimental Study of the Propagation 
of Sound Over Ground," J. Acoust. Soc. Am., Vol. 31, 1959, pp. 724- 
733. 

10. P.H. Parkin and W.E. Scholes, "The Horizontal Propagation of Sound 
From a Jet Engine Close to Ground, at Radlett," J. Sound Vib., 
Vol. 1, 1964, pp. 1-13. 

11. P.H. Parkin and W.E. Scholes, "The Horizontal Propagation of Sound 
From a Jet Engine Close to Ground, at Hatfield," J. Sound Vib., 
Vol. 2, 1965, pp. 353-374. 

12. S.P. Pao and L.B. Evans, "Sound Attenuation Over Simulated Ground 
Cover," J, Acoust. Soc. Am. , Vol. 49, 1971, pp. 1069-1075. 

13. H.G. Jonasson, "Sound Reduction by Barriers on the Ground," J . 
Sound Vib., Vol. 22, 1972, pp. 113-126. 


90 


14. A. Cheng, Unpublished manuscript. Mechanical and Industrial Engin- 
eering Department, University of Utah, Salt Lake City, Utah. 

15. W.K. Van Moorhem, "An Investigation of Acoustic Propagation in a 
Thermally Inhomoheneous Atmosphere," UTEC ME82-015, Mechanical and 
Industrial Engineering Department, University of Utah, Salt Lake 
City, Utah. 

16. W.K. Van Moorhem, "Acoustic Propagation in a Thermally Stratified 
Atmosphere," "Fourth Seminannual Report to the National Aeronau- 
tic and Space Agency-Grant NAG-1-283," UTEC ME 84-075, Mechanical 
and Industrial Engineering Department, University of Utah, Salt 
Lake City, Utah. 

17. J. Mathews and R.L. Walker, Mathematical Methods of Physics, (W.A. 
Benjamin, Inc., New York), 1964. 

18. L.P. Smith, Methods for Scientists and Engineers, (Prentice-Hall, 
Inc., New York), 1953. 

19. D.A. Sachs and A. Silbiger, "Focusing and Refraction of Harmonic 
Sound and Transient Pulses in a Stratified Media," J. Acoust. 
Soc. Am. Vol . 49, 1971, pp. 824-840. 

20. D.W. While and M.A. Pedersen, "Evaluation of Shadow Zone Fields by 
Uniform Asymptotics and Complex Rays," J. Acoust. Soc. Am. Vol. 69, 
1981, pp. 1029-1058. 

21. C.I. Chessell, "Propagation of Noise Along a Finite Impedance 
Boundary," J. Acoust. Soc. Am .. Vol. 62, 1977, pp. 825-834. 

22. J.E. Piercy and R.F. Embleton, "Noise Testing of Vehicles-Acoustic 
Propagation Phenomena," S.A.E. Special Publication 456, 1980, pp. 
13-19. 

23. J.E. Piercy, R.F. Embleton, and L.C. Sutherland, "Review of Noise 
Propagation in the Atmosphere," J. Acoust. Soc. Am., Vol. 61, pp. 
1403-1418. 

24. S.P. Pao, A.R. Wenzel and P.B. Oncley, "Prediction of Ground Effects 
on Aircraft Noise," NASA Technical Paper 1104. 

25. W.K. Van Moorhem, personal communication. 


APPENDIX C 


The propagation of plane waves in a thermally stratified 
atmosphere 
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The behavior of a plane wave reflecting from a finite impedance surface with a realistic 
atmospheric temperature profile is investigated. An approximate solution lias been implemented 
on a digital computer and this solution is presented graphically for a number of cases at low 
incidence angles. 

PACS numbers: 43.28.Kt, 43.28.Fb, 43.20.Bi 


INTRODUCTION 

The atmosphere is not isothermal under normal condi- 
tions, yet most acoustic propagation models make that as- 
sumption. Temperature inhomogeneity can occur as ran- 
dom fluctuations resulting from convection cells and 
turbulence or as slow diurnal variations in the atmosphere 
immediately adjacent to the ground surface. It is the latter 
case that this paper addresses. 

During daylight hours when sufficient insolation is 
present, the atmosphere is typically warmer near the ground 
surface than far above it, a lapse condition. Although these 
regions of strong temperature gradients extend, at most, 
only a few meters above the ground, most acoustic receivers 
are also in .his range of altitude. Ground mounted micro- 
phones, which are sometimes used to avoid ground reflec- 
tions, are unfortunately located so as to receive the maxi- 
mum effects of refraction due to temperature gradients. 

I. TEMPERATURE PROFILE 

Any study of the effects of thermal gradients on near- 
earth sound propagation requires the selection of a tempera- 
ture profile which is both realistic and mathematically trac- 
table. Geiger 1 gives some experimental results which 
indicate that the variation of temperature with altitude be- 
haves as the logarithm cf the height above the ground. Ten- 
ekes and Lumley 2 give an argument based on turbulence 
theory that leads to the same conclusion. 

Logarithmic profiles are difficult to deal with math- 
ematically. Most previous work (e.g., Pekeris 3 ) on the prob- 
lem of acoustic propagation in a stratified atmosphere have 
considered the profile 

r=7’ 0 /(l+az), (1) 

which yields an exact solution to the resulting differential 
equation. The vertical temperature profile used here 



is similar to Eq. (1) in its behavior near the ground (z = 0) in 
that it allows a steep gradient, but, more realistically, asymp- 
totically approaches a constant value high above the ground. 
Using data obtained by Willshire 4 and by Butterworth 5 and 
fitting Eq. (2) to this data, typical values for a ranged from 
1.2 to 5.9 m' 1 and for A T /T m from 0.003.. ,o 0.0342 for 
lapse conditions. Similar values of a and magnitude of 


I \ 


A T /T m were obtained for inversion conditions. In other 
cases, particularly where e lapse condition is being replaced 
by an inversion or vice versa, the fit can be very poor. 


II. ANALYSIS 


Consideration of the problem of propagation of plane 
waves in a thermally inhomogeneous atmosphere yields con- 
siderable insight into the more complex problem of propaga- 
tion of waves from a point source. The problem is also inter- 
esting in its own right ts a model of the propagation of waves 
from a point source when the source is very far from the 
refracting region. 

The problem considered is governed by the equation 


1 ?P 


+ V 2 /> = 0, 


c 2 ^) dt 2 
with the sound speed c[z) given by 

At the surface (2 = 0) the boundary condition 
— Zw — P, 


(3) 

(4) 

(5) 


where w is the vertical component of the acoustic fluid veloc- 
ity and Z the acoustical impedance of the surface, is applied. 
The solution is assumed to behave as 


P= exp 


i{o)t —x cos 6 + — z sin 6 j J 


-I- R (0 )exp 


iiat — x cos 0 — z sin 6 ) 

L \ c„ c„ )\ 


( 6 ) 


as z —* 00 . This condition implies the solution behaves as 
plane waves in a homogeneous atmosphere outside the re- 
gion of strong gradients. The subscript 00 , indicating evalua- 
tion as z — * 00 , is omitted below to simplify notation. The 
subscript zero will be used to indicate evaluation at the 
ground surface. 

Taking 

P = F (z)exp| 1 [cot — (a)/c]x cos 6 ] J (7) 

and substituting into Eq. (3) yields 


d 2 F a) 2 / 1 + az 

~dF + "7 \ 1 +az + AT/T 


— cos 2 6 


y 


=0, 


( 8 ) 


a second-order differential equation with a turning point at 
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( 22 ) 


z Tr =[\/a)[\AT/T) cot 2 0-1]. (9) 

Nayfeh 6 gives a method of obtaining an approximate solu- 
tion of this equation which is valid when y = a>/(ac)> 1 . This 
requirement is the same as that found for Eq. (3) to be valid, 
so no new approximations are required. The approximate 
solution can be expressed as 


F = 


_ *, 

[*W /J 




where 


[ff'Ml 


’/W) 


1/2 


( 10 ) 


rfc)= (Jr) 273 ^). 


( 11 ) 


ft 


*.= 


Ae* w/ * 


e* w/ *-iR 


0<cosJ<[l + AT/T)- 112 , 
cos0>(\ + AT/T)~ 1 ' 2 
and z^Zj-pf 

( 12 ) 

cos 0>(1 + AT/T)~ m 
and 0<z<Zn>, 


'A-R, 
Ae-™ 
e* r/6> — lR’ 


*2 = 


ARe* r/6> 
e* ,/6> — iR ’ 


O<cos0<(l + AT/T)~ V2 , 

cos 0X1 +dr/7’)- ,/2 

and z>z 7T , (13) 

cos 0X1 + AT /T)~ xl2 
and 0<z<. rp , 


A = yFv73sin' ,2 [0)r' ,<, a ,,2 e* iw/l2 \ (14) 


R= - 


^iM°)) + '# i (v(Q)) 
r/r 2 (i7(0)) -(- '#2 (t/( 0)) 
i Z g'(0) 


r = ay — 

2p 0 c o g'(0) 

tt> = (Z/r o rotf n g'(0)(l) 2 ", 


(15) 

(16) 

(17) 



-±4U L ta (I±ft 

2 T sin 0 V 1 — <f>) 

/ sin 2 0(1 -t -az + AT /T)-AT/T V n 
V sin 2 0(1 + az + AT/T) / ’ 


(18) 

(19) 


*,(£)« (!) , ^ ,/2 ^%(U J/J ). (20) 

and 

( 21 ) 

Values of the modified one-third order Hankel func- 
tions, A, and A 2 , have been tabulated and their properties 
discussed. 7 For the numuical analysis the functions A, and 
A 2 have been expressed in terms of Airy functions, Ai and Bi, 
where 


and 

A 2 (£) = **[Ai( — £)-MBi( — £)). (23) 

where k = ( 12) ,/6 e ~ " /4 and k • is the complex conjugate of 
k. 

The asymptotic expressions for the modified Hankel 
functions are available, 7 and if the parameter y is sufficiently 
large the usual result for the reflection of plane waves from a 
plane surface in a homogeneous atmosphere is obtained. 

The form and behavior of the solution depends upon 
whether the waves which are initially traveling downward, 
toward the ground, are reflected upward from the ground or 
if refraction turns the waves upward before the ground is 
reached. If cos 0<[(1 + dT/T)] -172 the turning point is be- 
low the ground surface and reflection from the ground sur- 
face occurs. In the opposite case the turning point is above 
the ground surface and the wave is refracted upward before 
the ground is reached. 

Based on meteorological data whered T /T ranged up to 
0.0342, rays with incidence angles less than about 10* could 
be refracted upwards before reflection from the ground 
could occur. Although these angles are quite small, they are 
within the range of common “look angles” for aircraft noise 
testing and refractive effects can be expected to play a signifi- 
cant role in this case. 

If the waves are refracted into upward propagation at 
the turning point, the solution depends u*x>n whether the 
receiver is above or below the turning point. Above the turn- 
ing point the solution describes an interference pattern simi- 
lar to that in a uniform atmosphere. Below the turning point 
exponential growth or decay occurs for large values of rj. 

The first term in Eq. (10) decays as the receiver moves 
downward, away from the turning point, if V — 1 is taken to 
equal — i and, similarly, the acoustic field below the turning 
point is reduced as 0 goes to zero and the turning point 
moves upward, if ln( — 1) = — in is chosen. The second 
term in Eq. ( 10) also decays with decreasing 0, but grows as z 
decreases from the turning point value. This behavior occurs 
since the second term represents the reflection from the 
ground surface of the exponentially decreasing first term. 

III. RESULTS 

Numerical evaluation of the solution for the propaga- 
tion of plane waves in a temperature gradient has been per- 
formed on a digital computer. Figure 1 is a typical example 
and is intended as the referent as parameters are varied in 
other drawings. Surface impedance values given by Ches- 
seU’s® relationship for grass Lave been used, with specific 
flow resistance chosen as 300 cgs units. A value of y = 10 
corresponds to a frequency of approximately 1000 Hz with 
a = 2 ra - 1 . Although this value of y and the additional value 
discussed below arc not extremely large they are given as 
examples and experiments will be required to demonstrate 
the applicability of the theory of these values. Figures 1 to 3 
show plots of height above the ground versus sound pressure 
level in decibels. The reference level is assumed equal to the 
amplituc e of the incoming wave at an infinite height. Thus 
the plots show the difference between the sound pressure 
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level of the incident wave in an isothermal atmosphere which 
is neither reflected nor refracted and the levels resulting 
from the refraction/reflection process in a nonisothermai 
atmosphere. 

For the conditions used to calculate Fig. 1(a), y = 10, 
AT/T = 0.2, 6 = 5*, and a = 2 m“ 1 the turning point and 
caustic is located at a height of about 0.8 m. Note that the 
maximum sound pressure level occurs somewhat above this. 
Above the turning point an interference pattern occurs with 
a slowly decreasing sound level in successively higher maxi- 
ma. Since the ratio of amplitudes of the incident and refract- 
ed waves is near unity, values of the sound levels above ap- 
proximately 6 dB must be due to increased amplitude of the 
waves beyond the values occurring far above the ground. 
The similar amplitude of the two waves also leads to the 
strong cancellation observed in the interference pattern. Be- 
low the turning point a rapid decay occurs as the shadow 
region is entered. 

Figure 1(b) is for the same conditions but with an isoth- 
ermal atmosphere, A T /T = 0. In this case the reflected wave 
is somewhat weaker than the incident wave and this is appar- 
ent from the levels reached in the interference maxima and 
minima. Comparison of Fig. 1(a) and (b) clearly points out 
the increased sound levels and variable spacing of the inter- 
ference pattern in the thermally inhomogeneous case. 

Figure 2(a) shows a similar case but with the incidence 
angle increased to 10*. In this case the turning point would be 
below the ground surface and reflection at the ground oc- 
curs. The reflected wave is weaker than the incident wave in 



FIG. I . Height versus sound pressure level for: (a) y — 10 ,AT/T„ = 0.02, 
0=5*, a = 2m ';(b)r = \0,AT/T m = 0, 0 = 5*. a = 2 m'. 



FIG. 2. Height versus sound pressure level for: (a) f = 10 ,AT /T m — 0.02, 
0- 10", a = 2 m-' ; (b)r= 10. AT/T m = 0.02. 0 = 2*. u = 2 m -1 . 



FIG. 3 Height versus sound pressure level for: (a) y = 10, A T/T m = 0.01, 
0=5*, a = 2 tn-';(b)r= 10.AT/T. = 0.02, 0 = 5*. a = 1 m"\ 
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this case and this is indicated by the levels occumng in the 
interference pattern. Figure 2(b) again shows a similar case 
but with Q — 2*. In this case the turning point is at r height cf 
about 7.7 m above ground and the rapid dec*, below the 
turning point is more appare; t than in Fig. 1(a). 

In Fig. 3 U)AT /T has been decreased to 0.01 yielding a 
turning poin at about 0. 1 5 m. Comparison to Fig / "i! shows 
that this change has reduced the envelope of the sound level 
maxima, and that the interference minima are at a higher 
sound level, both due to a slightly weaker rejected wave 
than in Fig. 1(a), but stronger than the isotherms' case shown 
in Fig. 1(b). Figure 3(b) has the parameter a decreased to 1 
m -1 spreading the region of significant temperature vari- 
ation upward and decreasing the temperature gradient. The 
parameter A is held constant at ten but since a was decreased 
this corresponds to a decreased frequency as is apparent 
from the interference pattern. The turning poin' for this case 
is at a height of about 1.6 m. 

IV. CONCLUSIONS 

The situation investigated hei *, pine w. ves approach- 
ing a finite impedance boundary at an arbitrary angle 
through a lapse temperature gradient is se-.-n .o demonstrate 
the major phenomena occurring at ^ it pe distance from a 
point source. The interference pattern is seen to be signifi- 
cantly distorted from the iso then' . ' situation by an inhomo- 


geneous layer about 0.3 m thick (a = 2 m ') and with a 
temperature increase of about 6 *C near the ground {AT / 

T = 0.02) . 
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